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Abstract 

Aircraft  with  stores  can  exhibit  aeroelastic  limit-cycle  oscillation  in  the  tran¬ 
sonic  regime.  Limit-cycle  oscillation  is  an  aeroelastic  phenomenon  characterized  by 
limited  amplitude,  self-sustaining  oscillations  produced  by  fluid-structure  interac¬ 
tions.  In  order  to  study  this  phenomenon  for  practical  configurations,  a  first-order 
accurate  code  was  developed  to  interface  a  modal  structural  model  with  a  com¬ 
mercial,  parallel,  Euler /Navier- Stokes  fluid  solver  with  a  deforming  grid  capability. 
Initial  testing  of  this  code  was  completed  on  the  Goland"1"  and  AGARD  445.6  wings. 

Limit-cycle  oscillation  was  simulated  for  the  Goland"1"  rectangular  wing  and  was 
successfully  compared  to  previously  reported  findings.  It  was  found  that  the  aero¬ 
dynamic  nonlinearity  responsible  for  limit-cycle  oscillation  in  the  Goland"1"  wing  was 
the  periodic  appearance/disappearance  of  shocks.  The  Goland"1"  structural  model 
was  such  that  in  the  transonic  flutter  clip  region,  the  primary  bending  and  twisting 
modes  were  in  phase  and  coupled  to  produce  a  single-degree-of-freedom,  torsional 
flutter  mode  about  a  point  located  ahead  of  the  leading  edge  of  the  wing.  It  was 
determined  that  the  combination  of  strong  trailing-edge  and  lambda  shocks  which 
periodically  appear/disappear,  limited  the  energy  flow  into  the  structure.  This  mech¬ 
anism  quenched  the  growth  of  the  flutter,  resulting  in  a  steady  limit-cycle  oscillation. 

Under- wing  and  tip  stores  were  added  to  the  Goland"1"  wing  to  determine  how 
they  affected  limit-cycle  oscillation.  It  was  found  that  aerodynamic  store  shapes 
affect  limit-cycle  oscillation  in  two  offsetting  ways:  by  interfering  with  the  flow  field 
on  the  wing  surface,  and  by  transferring  additional  store  forces  into  the  structure. 
The  under-wing  stores  interfere  with  the  airflow  on  the  lower  surface  of  the  wing 
which  decreases  limit-cycle  oscillation  amplitudes,  whereas,  under-wing  and  tip  store 
forces  transferred  into  the  wing  structure  directly  increase  limit-cycle  oscillation  am¬ 
plitudes. 
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I.  Introduction 

Aeroelasticity  is  the  interaction  between  an  clastic  structure  in  an  air  stream  and 
the  resulting  aerodynamic  force  [29].  The  performance  of  aircraft  is  often  limited  by 
adverse  aeroelastic  interactions  such  as  flutter.  Flutter  is  defined  as:  “A  dynamic 
instability  of  a  flight  vehicle  associated  with  the  interaction  of  aerodynamic,  elastic, 
and  inertial  forces”  [29].  Classical  flutter  is  characterized  by  catastrophic  diverging 
oscillations  [15]. 

Aircraft  operating  in  the  transonic  region  risk  encountering  limit-cycle  oscilla¬ 
tion  (LCO)  which  is  also  called  limit-cycle  flutter  and  limited  amplitude  flutter  [8]. 
Limit-cycle  oscillations  are  limited  amplitude,  self-sustaining  oscillations  produced 
by  fluid-structure  interactions.  LCO  results  in  undesirable  airframe  vibrations  that 
adversely  affect  a  pilot’s  ability  to  function  and  degrades  targeting  accuracy  for 
fighter  aircraft.  Aircraft  with  high  aspect  ratio  or  highly  swept  wings  can  experience 
LCO  [41],  For  typical  LCO,  the  amplitude  is  constant  for  a  given  flight  speed.  The 
amplitude  increases  when  the  flight  speed  is  increased  and  a  new  constant  amplitude 
is  obtained  when  the  flight  speed  is  again  fixed  [8].  If  the  flight  speed  is  decreased 
but  is  still  above  the  flutter  point,  a  new  constant,  lower  amplitude  is  obtained.  If 
the  flight  speed  is  decreased  below  the  flutter  point,  the  oscillations  damp  out.  For 
nontypical  LCO,  the  amplitude  does  not  increase  as  the  flight  speed  increases  [8]. 

LCO  is  related  to  flutter,  which  occurs  when  the  dynamic  pressure  is  increased 
to  a  critical  value  above  which  the  system  becomes  dynamically  unstable.  When  a 
non-linear  mechanism  that  counters  the  amplitude  growth  is  present  in  this  unstable 
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region,  LCO  can  result.  The  nonlinearity  that  sustains  LCO  can  be  in  the  struc¬ 
ture,  the  aerodynamic  flow,  or  in  both.  The  case  of  interest  for  this  study  was  the 
nonlinearity  in  the  aerodynamic  flow.  Sources  of  aerodynamic  nonlinearity  include 
separated  flow  and  shock  motion  in  transonic  flow  [17]. 

Advances  in  nonlinear  modelling  and  computer  hardware  have  led  to  nonlin¬ 
ear  aeroelastic  predictions  for  reasonably  complex  configurations  [22,  23,  24,  42,  44], 
However,  there  is  continued  interest  in  LCO  motivated  by  the  need  to  further  under¬ 
stand  the  physics  and  mechanisms  involved  in  LCO  so  that  better  LCO  predictive 
tools  can  be  developed  [12,  15,  52],  This  research  examines  the  aerodynamic  nonlin¬ 
earities  for  a  clean  wing  and  the  non-linear  aerodynamic  effects  induced  by  under¬ 
wing  stores.  There  have  been  several  efforts  where  LCO  driven  by  aerodynamic 
nonlinearities  were  studied  [35,  42,  44],  This  research  adds  to  this  knowledge  base 
by  providing  a  detailed  analysis  of  the  shock  motion  during  LCO  and  by  determining 
the  role  of  store  aerodynamics  in  LCO. 

Tij deman  and  Seebass  [56]  characterized  the  periodic  motion  of  shock  waves  on 
oscillating  airfoils  into  three  different  types.  In  type-A  shock  motion,  the  shock  wave 
moves  fore  and  aft  sinusoidally.  The  shock  strength  may  vary  as  the  shock  moves, 
but  the  shock  is  always  present  [56].  In  type-B  shock  motion,  the  shock  moves 
sinusoidally  as  in  type-A,  but  the  shock  disappears  during  a  part  of  its  backward 
motion  [56].  Type-C  shock  motion  occurs  when  a  periodic  shock  wave  leaves  the 
airfoil  and  continues  upstream  as  a  weak,  free  shock- wave  [56].  Bendiksen  [5]  has 
shown  that  if  a  transonic  flutter  mode  that  appears  to  be  a  singlc-degree-of-freedom, 
torsional  motion,  ahead  of  the  leading  edge,  has  a  phase  angle  between  bending  and 
torsion  less  than  approximately  10  degrees,  and  has  a  transition  from  a  type-A  to 
type-B  shock  motion,  it  is  likely  that  LCO  will  result.  This  type-B  shock  motion 
limits  the  flow  of  energy  from  the  fluid  to  the  structure,  resulting  in  LCO  [5]. 

Under-wing  and  tip  stores  can  affect  LCO  by  changing  the  inertial  properties 
of  the  wing,  and  by  changing  the  aerodynamics  of  the  wing.  If  a  store  mass  is 
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added  forward  of  the  elastic  axis  (EA),  the  flutter  speed  generally  increases,  which 
produces  a  stabilizing  effect  [32],  If  a  store  mass  is  added  aft  of  the  elastic  axis,  it 
has  a  destabilizing  effect  [32],  Exceptions  to  this  can  occur  in  the  transonic  flutter 
dip  region  [6,  7].  Therefore,  if  there  is  typical  LCO  at  a  given  airspeed  and  a  mass  is 
added  ahead  of  the  elastic  axis,  and  it  raises  the  flutter  speed,  the  LCO  amplitude 
will  decrease  because  the  given  fixed  airspeed  becomes  closer  to  the  flutter  speed. 
If  the  mass  is  moved  aft  of  the  elastic  axis  and  it  lowers  the  flutter  speed,  the  LCO 
amplitudes  will  increase. 

The  role  of  store  aerodynamics  in  LCO  is  less  well  understood.  One  can  hy¬ 
pothesize  two  primary  mechanisms  by  which  store  aerodynamics  affect  LCO:  the 
stores  could  change  or  interfere  with  the  airflow  on  the  wing,  thereby  changing  how 
the  shock  motion  is  limiting  the  energy  transferred  into  the  structure;  and  the  store 
aerodynamic  carriage  loads  that  are  transferred  into  the  structure  could  sufficiently 
change  the  total  forces  experienced  by  the  wing,  thereby  changing  the  LCO  response. 
The  first  mechanism  deals  with  changes  to  the  wing  aerodynamics  because  of  the 
presence  of  the  store.  The  second  mechanism  is  independent  of  any  changes  to  the 
wing  aerodynamics  and  any  changes  in  the  LCO  are  because  of  aerodynamic  forces 
experienced  by  the  store.  The  most  likely  method  that  store  aerodynamics  affect 
LCO  is  through  a  combination  of  these  two. 

1.1  Overview  of  Computational  Aeroelasticity 

Several  computational  methods  have  been  used  to  model  fluid-structure  inter¬ 
actions  such  as  flutter  and  LCO.  The  simplest  and  fastest  method  is  to  use  a  finite 
element  formulation  for  the  structure  and  a  linear  aerodynamic  model  such  as  the 
doublet-lattice  method  for  the  fluid  [31].  Finite  element  models  are  very  fast  and 
accurately  predict  instability  frequencies,  therefore,  they  are  often  used  as  a  first 
step  in  aeroelastic  analysis  [6,  7,  12,  13,  15,  51].  However,  in  the  transonic  region 
where  the  aerodynamic  flow  is  non-linear,  these  methods  do  not  accurately  predict 
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flutter  onset  velocity  [12].  Accuracy  can  be  improved  when  using  a  finite  element 
formulation  for  the  structure  by  employing  an  unsteady-aerodynamics  computation 
procedure  such  as  an  aerodynamic-influence-coefficient  matrix.  This  procedure  cal¬ 
culates  the  unsteady  aerodynamics  once  and  stores  the  results  in  tables  that  are 
reused  for  repetitive  aeroelastic  computations  [10]. 

Greater  accuracy  can  be  obtained  by  using  a  more  complex  aerodynamic 
model.  A  computationally  efficient  aerodynamic  model  that  is  frequently  used 
solves  the  three-dimensional  transonic  small  disturbance  potential  flow  equations. 
This  method  has  been  used  in  research  codes  [34]  and  in  the  NASA  Langley  de¬ 
veloped  CAPTSDv  [3,  6,  7,  14,  31,  51].  Mildly  separated  flows  have  been  studied 
with  CAPTSDv  with  the  addition  of  an  inverse  integral  boundary  layer  model  [31]. 
These  models  are  more  accurate  than  linear  aerodynamic  models  but  they  do  not 
completely  capture  the  intricacy  of  the  flow  due  to  shock-boundary  layer  interaction 
or  large  angles  of  attack. 

The  latest  aeroelastic  codes  solve  the  Euler  equations  for  inviscid  flows  or  the 
Navier-Stokes  equations  for  viscous  flows  for  even  greater  accuracy,  but  at  a  cost  of 
greater  computational  time.  The  most  widely  available  of  these  codes  are  ENS3DAE 
and  CFL3DAE.  ENS3DAE  was  developed  by  the  Lockheed-Georgia  Company.  It 
uses  a  central  finite-difference  scheme  and  structured  grids.  It  has  been  used  to  solve 
a  large  number  of  aeroelastic  problems  [6,  7,  31,  51].  CFL3DAE  was  developed  by 
NASA  Langley.  It  uses  an  upwind  finite-volume  scheme  and  structured  grids.  It  has 
also  been  used  to  solve  a  large  number  of  aeroelastic  problems  [3,  6,  7,  31,  51].  In 
addition  to  these  widely  distributed  codes,  the  majority  of  the  latest  research  codes 
use  a  finite-volume  scheme  to  solve  the  Euler/Navier- Stokes  equations.  Most  of  these 
codes  are  based  on  structured  grids  [26,  27,  37,  42,  43,  44,  45,  46,  53]  while  a  few 
recent  codes  use  unstructured  grids  [20,  21,  22,  23,  33,  35,  47,  48]. 

In  an  effort  to  attain  the  greater  accuracy  obtained  from  using  Euler/Navier- 
Stokes  solvers,  without  the  computational  expense,  there  is  significant  research  in 
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the  area  of  reduced  order  modelling  (ROM).  In  these  methods,  the  dominant  spatial 
modes  of  the  flow  held  are  determined  and  used  to  represent  the  flow  [16] .  The  most 
common  ROM  techniques  include  harmonic  balance  [54,  55]  and  proper  orthogonal 
decomposition  [16].  ROM  can  reduce  the  computational  expense  by  several  orders  of 
magnitude  but  care  must  be  taken  that  the  flow  held  is  not  over  simplified,  thereby, 
excluding  or  changing  the  nonlinearity  being  studied. 

Computational  methods  have  successfully  been  used  to  predict  the  flutter  point 
of  several  aircraft  wings.  The  most  common  case  studied  [9,  27,  33,  34,  35,  37,  38, 
40,  45,  47,  48,  49]  is  the  AGARD  445.6  weakened  wing,  which  has  been  presented 
as  a  standard  aeroelastic  test  case  [57].  Experimental  flutter  data  for  this  wing  is 
available  at  several  points  between  Mach  0.499  and  Mach  1.141  [57].  This  test  case 
has  often  been  used  to  validate  aeroelastic  programs.  With  various  computational 
aeroelasticity  methods,  good  agreement  with  experimental  data  was  achieved  for 
subsonic  results  but  all  the  computational  methods  tend  to  over  predict  the  flutter 
point  for  supersonic  results  [9,  27,  33,  34,  35,  37,  38,  40,  45,  49].  It  has  been  found 
that  inclusion  of  turbulence  and  viscous  effects  slightly  improved  the  flutter  point 
predictions  of  the  AGARD  445.6  wing  [27]. 

Computational  methods  have  also  been  successfully  used  to  predict  the  flutter 
point  of  fighter  wings.  Computational  results  compared  well  with  flight  test  results 
for  an  F-16  in  clean  configuration  [22,  23,  43]  and  with  stores  [10].  It  has  been  noted 
that  the  level  acceleration  achievable  by  an  F-16  does  not  significantly  impact  its 
aeroelastic  parameters.  However,  increasing  the  wing  loading  by  increasing  the  angle- 
of-attack  increased  the  aeroelastic  torsional  frequency.  For  a  generic  delta  fighter 
wing  with  stores,  it  has  been  found  that  adding  or  changing  the  store  in  the  structure 
affected  the  flutter  point,  however  store  aerodynamics  only  affected  the  supersonic 
flutter  point,  not  the  subsonic  [34,  53].  It  has  also  been  found  that  computational 
aeroelasticity  methods  which  take  into  account  aerodynamic  nonlinearities  produce 
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a  much  more  conservative  flutter  boundary  than  linear  theory  in  high  transonic 
flow  [33,  35]. 

All  of  the  previously  listed  computational  aeroelasticity  methods  have  been 
used  to  study  wings  with  linear  structural  models  that  exhibit  transonic  limit-cycle 
oscillations.  For  the  Goland+  rectangular  wing  with  a  tip  store  [6,  7,  51],  researchers 
found  that  LCOs  occurred  at  speeds  less  than  that  predicted  using  linear  analysis 
and  at  speeds  lower  than  that  computed  for  a  clean  wing.  They  also  found  that 
modelling  the  aerodynamics  of  the  tip  store  did  not  affect  the  LCO.  Each  of  the 
computational  methods  they  tested  showed  LCO  responses  over  the  Mach  range 
0.91  through  0.94  but  with  different  onset  velocities  [6].  The  computed  frequencies 
were  qualitatively  similar  but  did  not  agree  in  magnitude  as  well  as  expected  [6]. 

Attempts  have  been  made  to  model  LCO  in  the  F-16  in  both  a  clean  con¬ 
figuration  and  with  stores.  Melville  discovered  a  limit-cycle  oscillation  case  for  an 
F-16  without  stores  at  10  degrees  angle-of-attack  [42,  44],  He  determined  that  the 
shock-induced,  trailing-edge  separation  was  a  critical  nonlinear  limiting  mechanism 
responsible  for  LCO.  Denegri  has  identified  F-16  store  configurations  that  lead  to 
LCO  during  flight  test  [13,  15].  These  configurations  have  been  used  to  computa¬ 
tionally  model  LCO  with  mixed  results  [14,  47,  48,  54,  55].  The  computed  LCO 
amplitudes  have  not  successfully  duplicated  those  seen  in  flight  test.  Research  on 
this  configuration  is  continuing  in  order  to  determine  what  fidelity  of  computational 
model  is  required  to  correctly  simulate  this  LCO. 

1.2  Scope  of  Research 

The  objective  of  this  research  is  to  determine  how  nonlinear  aerodynamics 
quench  the  energy  flow  from  the  fluid  into  the  structure  resulting  in  LCO.  ft  is  known 
that  shock  motion  is  the  nonlinearity  responsible  for  LCO  in  the  Goland+  wing,  but 
it  is  not  understood  how  this  shock  motion  quenches  the  flutter  growth  resulting 
in  LCO.  It  is  also  known  that  adding  under- wing  stores  to  the  structure  affects  the 
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flutter  onset  speed.  However,  it  is  not  understood  how  the  store  aerodynamics  affect 
the  LCO  response  nor  is  it  understood  how  they  affect  the  quenching  mechanism 
responsible  for  LCO.  The  scope  of  this  research  is  defined  by  the  following  thesis 
statement  and  research  approach. 

1.2.1  Thesis  Statement.  The  primary  bending  and  torsional  modes  that  are 
known  to  couple  in  the  Goland+  wing,  result  in  periodic  shocks  that  are  responsible  for 
inhibiting  divergent  flutter  by  decreasing  the  restoring  forces  on  the  wing,  resulting 
in  a  balancing  of  the  inertial  forces  and  a  stable  limit-cycle  oscillation.  Further, 
store  aerodynamics  do  not  cause  limit-cycle  oscillation  but  can  decrease  or  increase 
the  amplitude  of  the  limit-cycle  oscillation  depending  on  the  size  and  location  of  the 
store. 


1.2.2  Research  Approach.  In  order  to  address  the  thesis  statement,  a  com¬ 
putational  aeroelasticity  code  was  developed  around  a  commercial  computational 
fluid  dynamics  (CFD)  code  and  a  linear  modal  structural  model.  It  was  required 
that  the  CFD  code  have  the  capability  to  interact  with  user  written  programs.  This 
interaction  allowed  a  structural  model  to  be  linked  to  the  CFD  code.  The  CFD  code 
also  had  to  incorporate  a  deforming  mesh  capability  in  order  to  simulate  the  flow 
around  a  moving  wing.  Finally,  the  CFD  code  had  to  use  unstructured  grid  method¬ 
ology.  Unstructured  grids  were  required  in  order  to  speed  up  the  grid  generation 
time  since  multiple  configurations  were  studied.  In  addition  to  these  requirements, 
it  was  desired  that  the  CFD  code  solve  the  Navier-Stokes  equations  and  have  multi¬ 
ple  turbulent  models  for  greater  accuracy.  FLUENT  6.1  by  FLLTENT  Incorporated 
was  chosen  because  it  met  these  requirements.  Once  the  CFD  code  was  chosen,  it 
was  linked  to  a  modal  structural  model  in  order  to  create  an  aeroelastic  program. 

Once  the  aeroelastic  program  was  developed  and  validated,  it  was  used  to 
simulate  transonic  LCO.  External  aerodynamic  stores  with  pylons  were  then  added 
in  order  to  analyze  how  the  LCO  response  changed.  The  stores  were  not  added  to 
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the  structural  model,  only  to  the  aerodynamic  model.  Modal  response,  force  and 
moment  data,  velocity  vectors,  pressure  contour  maps,  Mach  contour  maps,  density 
contour  maps,  power  spectral  density  plots,  work  on  the  structure,  and  wing  tip 
response  were  studied  in  order  to  quantify  how  the  store  aerodynamics  affected  the 
LCO  and  to  determine  what  aerodynamic  nonlinearity  was  providing  the  quenching 
mechanism  responsible  for  the  LCO. 

1.3  Document  Organization 

This  document  is  organized  into  five  chapters.  Chapter  I  introduces  the  prob¬ 
lem,  discusses  prior  work,  states  the  thesis  of  the  research,  and  lays  out  the  organi¬ 
zation  of  the  document. 

Chapter  II  provides  an  overview  of  FLLIENT  6.1.  It  also  describes  the  aeroe- 
lastic  program  that  was  written  to  interact  with  FLLIENT  in  order  that  it  could  be 
used  as  an  aeroelastic  analysis  tool. 

Chapter  III  describes  the  validation  of  the  aeroelastic  program.  It  covers  the 
prescribed  motion,  Goland"1"  wing,  and  AGARD  445.6  wing  test  cases  used  in  vali¬ 
dating  this  aeroelastic  analysis  method. 

Chapter  IV  discusses  the  results  obtained  using  the  aeroelastic  program.  The 
chapter  is  divided  into  three  sections.  The  first  section  discusses  the  Goland"1"  wing 
without  aerodynamic  stores  and  analyzes  the  quenching  mechanism  responsible  for 
the  LCO.  The  second  section  discusses  the  Goland"1"  wing  with  aerodynamic  under¬ 
wing  stores  and  their  effect  on  LCO.  The  last  section  discusses  the  Goland"1"  wing 
with  an  aerodynamic  tip  store  and  its  effect  on  LCO. 

Chapter  V  summarizes  the  conclusions  and  provides  recommendations  for  fu¬ 
ture  research. 

In  addition  to  these  five  chapters,  this  document  contains  supporting  appen¬ 
dices  that  provide  detailed  information  in  support  of  the  Eve  main  chapters. 
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II.  Aeroelastic  Program 


2.1  Aeroelastic  Analysis  Method 

One  goal  of  this  research  effort  was  to  extend  the  investigation  of  Dr  Beran, 
et  al.  [6,  7,  51]  by  evaluating  the  capability  of  a  commercial,  parallel,  Euler/Navier- 
Stokes  fluid  solver,  coupled  with  a  modal-analysis  structure  model,  to  accurately 
predict  aeroelastic  behavior  such  as  flutter  and  LCO.  There  is  not  currently  available 
a  dynamic-aeroelastic  analysis  tool  based  on  a  fully  unstructured  formulation.  This 
is  needed  to  reduce  the  grid-generation  time  in  order  to  investigate  LCO  for  complex 
wing/store/pylon  configurations.  Grid  generation  for  these  configurations  is  difficult 
for  structured  methods  and  usually  requires  overset  technology  [4,  39]. 

FLUENT  6.1  by  FLLTENT  Incorporated  was  chosen  to  fill  the  need  for  a  com¬ 
mercial  solver  that  could  handle  geometrically-complex  configurations  undergoing 
structural  deformation.  FLLIENT  6.1  is  an  unstructured  fluid  solver  with  a  deform¬ 
ing  grid  capability  that  can  be  controlled  through  a  user-written  subroutine  called  a 
user-defined  function  (UDF).  A  LIDF  was  written  that  extracted  aerodynamic  force 
data  from  the  fluid  solver  and  applied  it  to  a  modal  structure  model  in  order  to 
deform  the  grid.  By  coupling  FLUENT  with  a  modal  structure  model,  groundwork 
was  laid  to  analyze  aeroelastic  problems  for  complex  geometries. 

2.2  Overview  of  FLUENT  6.1 

FLUENT  is  a  commercially  available  computer  program  for  modelling  fluid 
flow  and  heat  transfer  in  complex  geometries.  FLLIENT  solves  flow  problems  using 
unstructured  meshes  that  can  be  generated  about  complex  geometries  with  relative 
ease  [2],  FLUENT  uses  a  control- volume-based  technique  which  consists  of  integrat¬ 
ing  the  governing  equations  on  a  control- volume  basis  yielding  discrete  equations  [2], 
FLLIENT  has  two  different  fluid  solvers:  a  segregated  solver  and  a  coupled  solver. 
The  segregated  solver  is  for  low  speed  or  incompressible  fluids.  For  compressible 
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fluids,  the  coupled  solver  should  be  used.  The  coupled  solver  solves  the  govern¬ 
ing  equations  of  continuity,  momentum,  energy,  and  species  transport  (if  required) 
simultaneously.  The  governing  equations  can  be  written  in  vector  form  as: 


—  j  WdV  +  j(F-G)-dA=  j  HdV 


where  W,  F,  and  G  are  defined  as 


(2.2) 


and  H  contains  source  terms  such  as  body  forces  and  energy  sources  [2], 

Because  the  governing  equations  are  non-linear  and  coupled,  several  iterations 
must  be  performed  to  obtain  convergence.  Each  iteration  consists  of  the  following 
steps  [2]: 


1.  Fluid  properties  are  updated,  based  on  the  current  solution.  (If  the  calculation 
has  just  begun,  the  fluid  properties  will  be  set  to  the  initial  conditions.) 

2.  The  continuity,  momentum,  and  energy  equations  are  solved  simultaneously. 

3.  Where  appropriate,  equations  for  scalars  such  as  turbulence  are  solved  using 
the  previously  updated  values  of  the  other  variables. 

4.  A  check  for  convergence  is  made. 

These  steps  are  continued  until  the  convergence  criteria  is  met  [2],  When 
performing  unsteady  simulations,  FLUENT  uses  an  implicit-time  formulation  (dual¬ 
time  stepping)  which  treats  each  physical  time  step  of  the  unsteady  problem  as 
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a  steady  state  problem,  which  is  run  to  convergence.  The  solver  can  be  run  with 
first  or  second-order  temporal  accuracy,  however,  the  deforming  mesh  algorithms  are 
only  supported  by  the  first-order  time-accurate  solver.  The  coupled  solver  employs 
a  first  or  second-order,  spatially-accurate,  upwind  method  to  solve  the  governing 
equations  [2], 

FLUENT  includes  three  mesh-motion  methods  that  can  be  used  to  update 
the  volume  mesh  in  the  deforming  regions  of  the  grid  based  on  motion  defined  at 
the  boundaries:  spring-based  smoothing,  dynamic  layering,  and  local  remeshing  [2]. 
In  the  spring-based  smoothing  method,  the  edges  between  any  two  mesh  nodes  are 
idealized  as  a  network  of  interconnected  springs.  The  initial  configuration  of  the 
edges  before  any  boundary  motion  constitutes  the  equilibrium  state  of  the  mesh  [2]. 
A  displacement  at  a  given  boundary  node  will  generate  a  force  proportional  to  the 
displacement  along  all  the  springs  connected  to  the  node.  The  second  method, 
dynamic  layering,  can  be  used  in  prismatic-mesh  zones  to  add  or  remove  layers  of  cells 
adjacent  to  a  moving  boundary  [2],  The  third  method  is  remeshing.  When  the  spring- 
based  smoothing  method  is  used,  the  boundary  displacement  may  become  large 
compared  to  the  local  cell  sizes.  This  can  cause  the  cell  quality  to  deteriorate  or  the 
cells  to  become  degenerate  [2],  This  can  result  in  negative  cell  volumes  or  convergence 
problems  when  the  solution  is  updated  to  the  next  time  step.  To  circumvent  this 
problem,  FLLTENT  agglomerates  cells  that  violate  user-specified  skewness  or  size 
criteria  and  locally  remeshes  the  agglomerated  cells  [2],  If  the  new  cells  satisfy  the 
skewness  and  the  size  criteria,  the  mesh  is  locally  updated  with  the  new  cells  [2], 

FLLTENT  provides  several  turbulent  models  for  modelling  turbulent  flows,  in¬ 
cluding  the  Spalart-Allmaras  model,  standard  k  —  e  model,  renormalization-group 
(RNG)  k  —  e  model,  realizable  k  —  e  model,  standard  k  —  u  model,  shear-stress  trans¬ 
port  (SST)  k  —  uj  model,  Reynolds  stress  model  (RSM),  and  large-eddy  simulation 
(LES)  model  [2], 
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Figure  2.1  Aeroelastic  Program  Flowchart 
2.3  Overview  of  Aeroelastic  Program 

The  aeroelastic  program  developed  around  FLUENT  is  depicted  in  Figure 
2.1.  FLUENT  solves  the  Euler/Navier-Stokes  equations  (equations  2.1  and  2.2) 
and  calculates  flow  variables  such  as  velocity,  pressure,  density,  temperature,  etc., 
for  each  cell  in  the  flowfield.  It  also  deforms  the  mesh  based  on  new  wing-node 
coordinates  when  they  are  provided.  The  remainder  of  the  aeroelastic  program  is  a 
UDF  (written  in  C)  that  is  called  by  FLUENT.  In  addition  to  the  aerodynamic  and 
structural  grids,  the  aeroelastic  program  has  five  main  inputs:  two  spline  matrices, 
a  mode  shape  matrix,  a  mass  matrix,  and  a  stiffness  matrix.  The  mode  shape 
matrix,  the  mass  matrix,  and  the  stiffness  matrix  are  determined  during  modal 
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analysis  preprocessing  using  software  such  as  MSC/NASTRAN.  The  output  from 
preprocessing  is  then  reformatted  for  input  into  the  aeroelastic  program. 

In  general,  the  grid  used  by  the  structural  model  and  the  aerodynamic  grid 
used  by  FLUENT  are  non-point  matching.  A  spline  matrix  [S]  is  used  to  interpolate 
displacement  data  from  the  structural  grid  [5S]  to  the  aerodynamic  grid  [<5a]  using 
the  equation 

[<y  =  [SIN-  (2.3) 

It  is  also  used  to  interpolate  force  data  from  the  aerodynamic  grid  [Fa]  to  the  struc¬ 
tural  grid  [Fs\  using  the  equation 


[Fs]  =  m^].  (2.4) 

Because  FLLIENT  provides  aerodynamic  forces  at  cell  centers  rather  than  at 
grid  nodes,  the  coupling  between  the  fluid  and  structural  models  is  accomplished  with 
two  splines.  The  first  spline  matrix  (Si)  is  created  to  spline  displacements  between 
the  structural  grid  points  and  the  aerodynamic  grid  points.  FLUENT  calculates 
pressures  at  the  cell  centers.  These  pressures  are  then  extrapolated  to  the  face  cen¬ 
troids  on  the  wing  and  used  with  the  wall-face  area  vectors  to  calculate  the  pressure 
forces.  These  pressure  forces  are  then  combined  with  the  viscous  shear  forces  to  get 
the  total  forces.  Since  the  total  forces  are  at  the  face  centroids  and  not  at  the  aero¬ 
dynamic  grid  points,  a  second  spline  matrix  (S'2),  is  created  to  spline  between  the 
structural  grid  points  and  the  wall-face  centroids  on  the  aerodynamic  grid.  Coordi¬ 
nates  for  the  structural  grid  nodes,  aerodynamic  grid  nodes,  and  aerodynamic-grid 
wall-face  centroids  are  used  during  preprocessing  to  create  these  two  spline  matrices. 

The  two  spline  techniques  chosen  for  this  application  were  the  thin- plate  spline  [1, 
18,  50]  and  the  infinite-plate  spline  [1,  28,  50].  Details  of  the  spline  algorithms  are 
presented  in  Appendix  A.  Both  the  thin-plate  and  infinite-plate  spline  are  based  on 
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the  undeformed  node  coordinates  of  the  two  computational  grids,  therefore,  they  are 
computed  and  stored  during  pre-processing. 

The  Newmark  algorithm  [30],  presented  in  Appendix  B,  was  used  to  discretize 
the  structural  model.  The  semi-discrete  equations  of  motion  can  be  written  as 

+  [<?]{«} +  [*]{«}  =  {F}  (2.5) 

where  [M\  is  the  mass  matrix,  [C]  is  the  damping  matrix,  [K]  is  the  stiffness  matrix, 
{F}  is  the  vector  of  applied  forces,  and  {«},  {h},  and  {ii}  are  vectors  of  modal  dis¬ 
placements,  modal  velocities,  and  modal  accelerations,  respectively  [30].  Assuming 
a  linear  acceleration  and  no  structural  damping,  the  Newmark  algorithm  is  used  to 
solve  equation  2.5  for  {u]n+ i  with  the  following  set  of  equations: 

f  A  t2  \  f  At2  \ 

f  [M]  +  —[K}  \  {u}n+ 1  =  {Fjn+i  -  [K\  f  Mn  +  A t{u}n  +  —  {u}n\  ,  (2.6) 

{h}n+i  =  {u}n  +  —  (({'h}n  +  {w}n+l)  ,  (2.7) 

and 

At2 

{u}n+ 1  —  {u}n  +  A t{u}n  H  —  (2{-u}n  +  {-uln+i)  .  (2-8) 

6 

FLUENT  uses  an  implicit,  dual-time  stepping  formulation  for  unsteady  flow 
simulations.  Because  of  this  dual-time  stepping  formulation,  and  the  lack  of  a  pro¬ 
vision  to  extract  wall-force  values  and  a  fractional,  physical,  time-step  size  during 
the  sub-iterations,  the  fluid  solver  and  structural  model  can  only  be  coupled  at  the 
end  of  each  physical  time-step.  Therefore,  the  structure  always  lags  the  fluid  solver 
by  one  time-step. 

The  aeroelastic  program  moves  an  attached  store  as  a  rigid  body  based  on 
how  the  pylon  attachment  points  move  on  the  wing  surface.  The  aerodynamic  forces 
and  moments  acting  on  the  store  are  calculated  and  transferred  into  the  structure 
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through  the  pylon,  in  order  that  any  aerodynamic  effects  on  the  store  are  felt  by 
the  structural  model.  These  forces  can  also  be  turned  off  or  ignored,  allowing  the 
aerodynamic  interference  on  the  bottom  of  the  wing  dne  to  the  store  to  be  isolated. 
In  order  to  determine  how  much  of  an  effect  the  store  forces  are  having  on  LCO  as 
compared  with  how  much  of  an  effect  the  store  interference  with  the  wing  is  having  on 
LCO,  the  aeroelastic  program  can  also  magnify  the  store  forces  that  are  transferred 
into  the  structure.  In  addition,  these  magnified  forces  can  be  reversed  in  direction  in 
order  to  further  isolate  their  effect.  The  aeroelastic  analysis  program  also  allows  for 
inertial  forces  to  be  added  to  the  store  aerodynamic  forces.  To  add  inertial  forces,  a 
mass  and  mass  location  are  entered  into  the  aeroelastic  program.  The  acceleration 
vector  of  the  mass  is  calculated  from  the  store  motion.  This  acceleration  vector  is 
multiplied  by  the  mass  to  obtain  an  inertial  force  vector  which  is  converted  to  forces 
and  moments  at  the  store  attachment  location.  These  forces  are  then  transferred 
into  the  original  structural  model. 

A  typical  aeroelastic  analysis  is  started  by  computing  an  initial,  steady-state 
solution  for  the  undeformed  wing.  This  solution  is  used  as  the  starting  point  for  the 
unsteady,  deforming-grid  computations.  At  the  start  of  the  unsteady  run,  the  forces 
at  the  wall-face  centroids  on  the  wing  are  calculated.  These  forces  are  then  splined 
using  equation  2.4  to  the  structural-grid  nodes.  New  deformed  structural-grid  co¬ 
ordinates  are  then  calculated  using  the  Newmark  method.  These  new  structural 
displacements  are  then  splined  using  equation  2.3  to  get  new  aerodynamic-grid  coor¬ 
dinates.  FLLIENT  then  deforms  the  mesh  based  on  the  new  coordinates,  increments 
the  time-step,  and  calculates  new  flow  variables.  This  process  repeats  until  a  speci¬ 
fied  flow  time  is  reached. 

The  aeroelastic  program  is  fully  parallelized  to  take  advantage  of  using  mul¬ 
tiple  processors.  The  fluid  solver  is  second-order  accurate  spatially,  and  first-order 
accurate  temporally  (when  using  the  deforming  mesh  algorithms).  The  structural 
model  is  second-order  accurate,  both  spatially  and  temporally. 
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III.  Validation  of  Aeroelastic  Program 

The  aeroelastic  program  described  in  chapter  II  was  implemented  as  FLUENT  user- 
defined  functions  in  the  C  programming  language.  The  components  of  the  aeroe¬ 
lastic  program  (splines,  structural  model,  etc.)  were  tested  individually  to  verify 
correctness.  The  thin-plate  spline  was  then  tested  to  obtain  an  estimate  of  the  error 
induced  in  the  system  when  transferring  data  from  one  grid  to  the  other.  The  full 
aeroelastic  program  was  then  tested  in  two  stages,  first  with  prescribed  structural 
motion,  then  through  full  aeroelastic  structural  motion  using  two  different  3-D  wing 
configurations:  the  AGARD  445.6  wing,  flutter  case,  and  the  Goland+  wing,  LCO 
case. 

The  aerodynamic  forces  were  calculated  by  the  fluid  solver  for  the  aerodynamic 
grid.  These  forces  were  then  transferred  to  the  structural  grid  through  the  use  of  a 
spline  matrix.  Any  errors  introduced  in  the  force  distribution,  by  the  spline  matrix, 
would  result  in  an  inaccurate  deformation  of  the  wing.  In  order  to  quantify  the 
errors  induced  when  transferring  forces  from  the  aerodynamic  grid  to  the  structural 
grid,  the  thin-plate  spline  was  tested  first  with  a  model  problem  and  then  with 
the  Goland+  LCO  case.  The  model  problem  was  a  flat  plate.  The  “aerodynamic” 
grid  was  created  by  evenly  dividing  the  top  and  bottom  surfaces  of  the  plate  into 
six  by  four  cells  with  a  node  at  the  center  of  each  cell.  The  “structural”  grid  was 
created  by  evenly  dividing  the  top  and  bottom  surfaces  of  the  plate  into  three  by 
two  cells  with  a  node  at  the  center  of  each  cell.  A  distributed  force  was  then  applied 
to  the  plate  and  the  exact  force  at  each  node  was  calculated  for  both  grids.  The 
calculated  force  on  the  aerodynamic  grid  was  then  splined  to  the  structural  grid  and 
these  values  were  compared  with  the  exact  values  that  were  calculated.  Errors  at 
each  structural  node  varied  between  3%  and  16%.  However,  when  the  forces  on  the 
aerodynamic  grid  and  splined  structural  grid  were  summed,  the  total  forces  were  the 
same.  When  moments  about  a  point  were  calculated  for  both  the  aerodynamic  grid 
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and  the  splined  structural  grid  and  then  summed,  the  total  moments  were  the  same 
for  both  grids.  The  model  problem  showed  that  the  spline  conserves  the  total  forces 
and  moments  on  a  global  level,  but  that  errors  exist  on  a  cell-by-cell  basis. 

To  characterize  the  error  for  a  more  realistic  geometry,  the  total  forces  and 
total  moments  for  both  the  aerodynamic  and  structural  grids  were  calculated  for 
the  Goland+  Wing  (discussed  in  section  3.2.2)  at  each  time-step  during  an  LCO 
simulation  run.  The  average  error  was  0.0003%  except  for  the  total  moment  about 
the  y-axis  which  had  an  average  error  of  0.0078%.  The  total  work  was  also  calculated 
for  each  grid.  Work  was  calculated  by  taking  the  dot  product  of  the  change  in  force 
at  each  node,  and  the  displacement  of  each  node.  The  results  were  summed  to 
provide  a  total  work.  The  work  was  compared  at  each  time-step  for  both  grids.  The 
difference  between  the  total  work  calculated  on  the  structural  grid,  and  the  total 
work  calculated  on  the  aerodynamic  grid  varied  from  —0.1  to  0.17  N  ■  m  with  an 
average  of  0.0004  N  -  m  (Figure  3.1),  or  0.0017%.  Overall,  there  is  a  very  small  error 
in  the  total  moment  about  the  y-axis  and  in  the  work  induced  by  the  thin-plate 
spline  using  these  grids  for  this  problem. 

3.1  Prescribed  Structural  Motion  Test  Case 

The  Goland+  wing,  described  in  section  3.2.2,  was  sinusoidally  pitched  ±0.5°  at 
3  Hz.  Coefficients  of  lift  and  moment  were  plotted  and  compared  to  Beran  et  al.  [6,  7]. 
The  grid  consisted  of  678,657  tetrahedral  cells.  Beran  et  al.  [6]  used  non-match-point 
flow  conditions.  FLUENT  requires  match-point  flow  conditions.  Therefore  it  was 
necessary  to  change  the  freestream  temperature  and  pressure  conditions  in  order  to 
match  a  given  Mach  number,  freestream  velocity,  and  freestream  density.  The  Mach 
number  was  set  at  0.92,  the  freestream  velocity  was  set  to  400  ft/sec  (121.92  m/s), 
and  the  freestream  density  was  set  to  0.0023771  slugs/ ft3  (1.225  kg/m3).  For  7  = 
1.4  and  R  =  1716.16  J^0fR  ( R  =  287.05  j^),  this  gave  a  freestream  temperature  of 
78.5245°i?  (43.6247  K)  and  a  freestream  pressure  of  320.9412  lb/ ft2  (15, 366.746  Pa). 
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Goland+  Clean  Wing 
Mach  0.92,  600  ft/sec 


Figure  3.1  Difference  in  Work  (N  ■  m) 

As  shown  when  plotting  the  coefficient  of  moment  (Cm)  versus  the  coefficient 
of  lift  (Cl)  (Figure  3.2),  the  results  from  FLUENT  compare  well  to  ENS3DAE  and 
CAPTSDv.  There  was  a  slight  asymmetry  in  the  FLLTENT-based  results  attributed 
to  grid  asymmetries  introduced  as  the  mesh  was  deformed.  Overall,  the  results 
indicated  that  the  aeroclastic  program  produced  similar  results  to  other  software 
used  for  aeroelastic  analysis. 

A  test  was  performed  to  determine  the  computational  impact  of  moving  the 
grid.  The  amount  of  time  needed  per  iteration  was  compared  for  a  steady  state 
calculation  versus  an  unsteady  calculation  that  called  the  moving  grid  routine  every 
iteration.  All  the  force  and  structural  model  calculations  were  performed,  but  the 
grid  coordinates  were  not  changed  so  the  flow  for  both  solutions  was  identical.  The 
test  was  performed  on  two,  four,  eight,  and  sixteen  Athlon  MP  2600+  processors 
with  2  GB  RAM  and  a  Myrinet  network.  As  shown  in  Figure  3.3,  the  aeroelastic 
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Figure  3.2  Cl  -  Cm  Phase  Portrait 

code  with  deforming  grid  capability  had  degraded  computational  efficiency  when 
more  than  eight  processors  were  used. 

3.2  Aeroelastic  Structural  Motion  Test  Cases 

3.2.1  AGARD  445-6  Wing.  The  first  3-D  wing  to  be  studied  with  the  full 
aeroelastic  program  was  the  AGARD  445.6  weakened-wing  model  [57].  This  wind 
tunnel  model  has  been  presented  as  a  standard  aeroelastic  case  [57].  Experimental 
flutter  data  is  available  for  this  wing  at  several  points  [57]  between  Mach  0.499  and 
Mach  1.141.  This  test  case  has  often  been  used  to  validate  aeroelastic  programs  [27, 
38,  40], 

This  wing  consists  of  a  NACA  65A004  airfoil  in  the  streamwise  direction  with 
a  panel  aspect  ratio  of  1.6525  and  a  panel  taper  ratio  of  0.6576.  The  span  is  2.5  ft 
(0.762  m),  the  root  chord  is  1.833  ft  (0.5587  m),  and  the  tip  chord  is  1.208  ft  (0.3682 
m).  The  quarter  chord  line  is  swept  45  degrees  (Figure  3.4). 
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Number  of  Processors 


Figure  3.3  Computational  Impact  of  Moving  Grid 


Figure  3.4  AGARD  445.6  Planform 
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Figure  3.5  AGARD  445.6  Aerodynamic  Grid 


An  ASTROS  model,  developed  by  Kolonay  [36],  was  obtained  from  Gord- 
nier  [27]  and  used  in  this  test.  This  structural  model  produced  the  first  fourteen 
natural  vibration  modes  for  the  AGARD  445.6  wing. 

An  unstructured  aerodynamic  grid  consisting  of  177,778  tetrahedral  cells  (Fig¬ 
ure  3.5)  was  created  using  Gridgen.  A  calculation  of  the  flutter  point  was  then 
made  using  the  aeroelastic  program.  A  dynamic  pressure  was  chosen  (table  3.1)  and 
the  solution  was  computed  for  four  oscillations.  If  the  oscillations  were  growing,  a 
lower  dynamic  pressure  was  chosen  and  the  solution  recomputed.  If  the  oscillations 
were  damped,  a  higher  dynamic  pressure  was  chosen.  This  was  continued  until  the 
flutter  point  was  bounded  and  a  dynamic  pressure  was  found  that  produced  neutral 
oscillations.  The  point  of  neutral  oscillations  was  the  flutter  point. 

For  this  problem,  an  initial  displacement  perturbation  of  0.01  was  given  to 
mode  1.  A  time  step  of  0.0003  was  also  used.  A  time-step  convergence  study  showed 
that  larger  time  steps  caused  the  oscillations  to  grow  while  smaller  time  steps  did 
not  change  the  solution.  For  the  AGARD  445.6  wing,  computational  flutter  results 
are  compared  by  calculating  the  ratio  of  the  computational  flutter  dynamic  pressure 
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Tabic  3.1  AGARD  445.6  Test  Conditions  for  Mach  0.96 


q/<le 

Velocity 

Temperature 

Pressure 

0.8 

893  ft/sec 

(272 

m/sec) 

360°i? 

(200  K) 

76.02  lb/ ft2 

(3639.67  Pa) 

1.0 

998  ft/sec 

(304 

m/sec) 

450°i? 

(250  K) 

95.02  lb/ ft2 

(4549.59  Pa) 

1.2 

1094  ft /sec 

(333 

m/sec) 

540°i? 

(300  K) 

114.02  lb/ ft2 

(5459.50  Pa) 

1.4 

1181  ft/sec 

(360 

m/sec) 

630°i? 

(350  K) 

133.03  lb/ ft2 

(6369.49  Pa) 

Table  3.2  Flutter  Point  Results,  Mach  0.96 


Method 

q/qe 

Gordnier  coarse  viscous 

1.15 

Gordnier  medium  viscous 

1.12 

Gordnier  fine  viscous 

1.05 

Gordnier  coarse  inviscid 

1.02 

Gordnier  medium  inviscid 

0.96 

Gordnier  fine  inviscid 

0.84 

Parker  inviscid 

0.86 

Wind-Tunnel  Experiment 

1.00 

Lee-Rausch  et  al. 

0.89 

Vermeersch  et  al. 

1.47 

(q)  over  the  experimental  dynamic  pressure  (qe).  Using  the  aeroelastic  program,  at 
Mach  0.96,  a  flutter  point  of  q/qe  =  0.86  and  a  flutter  frequency  of  12.65  Hz  was 
calculated.  This  compares  well  with  Gordnier’s  fine  inviscid  grid  results  as  shown  in 
table  3.2  [27], 

3.2.2  Goland+  Wing.  The  second  3-D  wing  to  be  studied  was  the  Goland+ 
wing.  The  Goland"1"  is  a  variant  of  the  heavy  Goland  wing  developed  as  a  transonic 
flutter  test  case  by  Eastep  and  Olsen  [19].  Based  on  the  original  Goland  wing  [25], 
the  heavy  Goland  wing  has  increased  mass  to  ensure  applicability  in  the  transonic 
regime.  The  Goland"1"  version  of  the  heavy  Goland  wing  was  modelled  with  a  box 
structure  beam  to  allow  a  variety  of  store  attachment  options  [7].  This  wing  was 
selected  as  a  test  case  because  data  for  ENS3DAE  and  CAPTSDv  exists  for  this 
configuration  [6,  7,  51]. 
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The  Goland+  wing  is  rectangular  and  cantilevered  from  an  infinite  midplane. 
The  wing  semi-span  is  20  ft  (6.096  m)  and  the  chord  (c)  is  6  ft  (1.8288  m).  The 
thickness,  r,  is  0.04  ft  (0.01219  m).  The  elastic  axis  is  located  2  ft  (0.6096  m)  from 
the  leading  edge.  The  airfoil  section  is  constant  over  the  spanwise  extent  of  the  wing 
and  is  a  symmetric,  parabolic-arc  airfoil  given  by 

*  =  ±2r(l- -)(-).  (3.1) 

c  c 

As  reported  by  Beran  et  al.  [7]  and  Snyder  et  al.  [51],  the  Goland+  wing-box 
structure  (Figure  3.6)  is  made  up  of  shear  elements  representing  three  spars  evenly 
spaced  2  ft  (0.6096  m)  apart,  and  11  ribs  evenly  spaced  2  ft  (0.6096  m)  apart.  Each 
rib  and  spar  is  4  in  (0.1016  m)  high.  Each  of  the  20  cells  created  by  the  ribs  and  spars 
are  capped  with  an  upper  and  lower  wing  skin  membrane  element.  Rods  are  added 
on  the  top  and  bottom  of  each  shear  element  and  at  every  spar/rib  intersection. 
These  137  rod  elements  represent  60  spar  caps,  44  rib  caps,  and  33  posts.  Every 
element  is  modelled  using  the  fictional  material  properties  shown  in  table  3.3  [51]. 
These  properties  were  selected  to  match  the  structural  dynamic  characteristics  of 
the  wing-box  model  to  those  of  the  beam  model  of  the  heavy  Goland  wing  [19]. 
The  elements  were  sized  to  minimize  the  differences  between  the  first  three  natural 
frequencies.  The  resulting  element  dimensions  are  shown  in  table  3.4  [51].  The  mass 
properties  are  modelled  by  placing  lumped  masses  with  no  rotational  inertia  at  each 
grid  point  [51].  The  internal  ribs  were  modelled  with  lumped  masses  of  1.9650  slugs 
(28.6768  kg)  at  each  leading  edge  spar,  3.9442  slugs  (57.5609  kg)  at  each  center  spar, 
and  5.3398  slugs  (77.9280  kg)  at  each  trailing  edge  spar.  The  masses  on  the  root 
and  tip  ribs  are  half  those  of  the  internal  ribs.  In  order  to  match  the  LCO  results  as 
presented  by  Snyder  et  al.  [51],  a  tip  store  was  added  to  the  structure.  The  tip  store 
is  10  ft  (3.048  m)  long  and  is  modelled  as  a  series  of  rigid  bar  elements.  The  bar  is 
positioned  6  in  (0.1524  m)  outboard  of  the  wing  tip  and  extends  3  ft  (0.9144  m)  in 
front  of  the  wing  leading  edge  and  1  ft  (0.3048  m)  behind  the  trailing  edge.  The  tip 
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Figure  3.6  Goland+  Structural  Model 


Table  3.3 

Goland+  Material  Properties 

Parameter 

Value 

Young’s  Modulus,  E 
Shear  Modulus,  G 
Structural  Density  ps 

1.4976xl09  slugs/ ft2  (2.3525xlOn  kg/m2) 
5.616xl08  slugs/ ft2  (8.822xl010  kg/m2) 
0.000001  slugs/ ft3  (0.0005154  kg/m3) 

store  has  a  mass  of  22.498  slugs  (328.3313  kg)  located  1.75  ft  (0.5334  m)  forward  of 
the  elastic  axis  and  a  rotational  inertia  of  50.3396  slugs-ft2  (68.2509  kg-rn2)  [51]. 

Natural  vibration  modes  for  the  Goland+  wing  with  a  tip  mass  were  obtained 
from  MSC/NASTRAN  [6,  51].  The  first  six  modes  were  used  in  the  aeroelastic 
analysis.  These  modes  are  plotted  in  Figure  3.7,  which  shows  the  top  of  the  right 
wing  with  the  leading  edge  to  the  upper  right  and  the  trailing  edge  to  the  lower  left. 

The  aeroelastic  code  was  tested  on  a  coarse,  unstructured,  tetrahedral  grid 
made  up  of  68,949  cells  with  518  cell  faces  on  the  wing  itself.  This  grid  was  pur¬ 
posefully  created  very  coarse  in  order  to  decrease  the  run  times  for  initial  testing 
of  the  aeroelastic  program.  The  program  was  run  at  multiple  Mach  numbers  and 
velocities  in  order  to  calculate  the  flutter  stability  boundary.  The  Goland+  sirnu- 
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Tabic  3.4  Goland+  Element  Dimensions 


Elements 

Dimensions 

Upper  and  lower  wing  skin  thickness 
Leading  and  trailing  edge  spar  thickness 
Center  spar  thickness 

Rib  thickness 

Post  area 

Leading  and  trailing  edge  spar  cap  area 
Center  spar  cap  area 

Rib  cap  area 

0.0155/f  (4.7244xl0~3m) 
0.0006/t  (1.8288xl0“4m) 
0.0889 ft  (2.7097x10 ~2m) 
0.0347 ft  (1.0577xl0~2m) 
0.0008/i2  (7.4322xl0"5m2) 
0.0416/t2  (3.8648xl0_3m2) 
0.1496/i2  (1.3898xl0_2m2) 
0.0422/t2  (3.9205xl0"3m2) 

(d)  Mode  4  (10.8 
Hz) 


(e)  Mode  5  (16.3  Hz) 


(f)  Mode  6  (22.8  Hz) 


Figure  3.7  Mode  Shapes  for  Goland+  Wing  with  Tip  Store 
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Figure  3.8  Goland+  Stability  Boundary 


lations  were  started  by  giving  mode  2  a  modal  displacement  perturbation  of  0.02. 
These  results  were  then  compared  to  the  results  presented  by  Snyder,  et  al.  [51]. 
Snyder,  et  al.  [51]  used  non-match-point  flow  conditions  for  their  analysis.  For  each 
test  case,  they  chose  a  Mach  number  and  velocity,  and  fixed  the  far-field  density  and 
temperature  values  at  standard-day,  sea-level  conditions.  In  order  to  match  their 
flow  parameters  in  FLUENT,  the  specific  gas  constant  R  and  the  far-field  pressure 
were  varied  for  different  Mach  number  and  velocity  settings.  The  specific  heat  ratio, 
7,  was  held  fixed  at  1.4.  The  computed  stability  boundary  is  shown  in  Figure  3.8. 
For  an  extremely  coarse  grid,  it  compares  well  with  the  results  presented  by  Snyder, 
et  al.  [51]. 

At  600  ft/sec  (182.88  m/sec)  and  Mach  0.92,  the  dynamic  pressure  is  beyond 
the  flutter  point  and  was  shown  by  Snyder,  et  al.  [51]  to  result  in  LCO.  At  these 
conditions  and  with  a  time  step  of  0.001  seconds,  the  aeroelastic  program  did  produce 
LCO  with  a  frequency  of  3.37  Hz.  A  time-step  convergence  study  showed  that  larger 
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Figure  3.9  Goland+  Coarse  Grid  Amplitudes 

time  steps  caused  the  oscillations  to  damp  while  smaller  time  steps  did  not  affect  the 
solution.  It  took  3  hours  and  17  minutes  to  produce  one  second  of  data  using  eight 

2.2  GHz  opteron  processors.  Mode  1  and  mode  2  amplitude  time  histories  as  well 
as  Cl  time  history  are  shown  in  Figure  3.9.  These  tests  showed  that  the  aeroelastic 
program  was  capable  of  modelling  LCO. 

3. 3  Validation  Summary 

An  aeroelastic  program  based  on  a  fully  unstructured-grid  formulation  capable 
of  simulating  flutter  and  LCO  was  successfully  developed  by  integrating  a  modal 
structural  model  with  FLUENT  6.1.  The  individual  parts  of  the  aeroelastic  code 
such  as  the  splines  and  structural  model  were  validated  to  ensure  they  produced  the 
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correct  results.  The  transfer  of  forces  from  the  aerodynamic  grid  to  the  structural 
grid  was  critical  for  the  accurate  calculation  of  wing  deformations.  The  thin-plate 
spline  was  therefore  tested  for  both  a  flat  plate  and  with  the  Goland+  wing  to  obtain 
an  estimate  of  the  errors  induced  by  the  spline.  It  was  found  that  there  was  an  error 
induced  but  on  average  it  was  less  than  0.0078%. 

Prescribed  motion  results  showed  that  the  aeroelastic  program  produced  lift 
and  moment  results  similar  to  those  produced  by  other  aeroelastic  analysis  codes. 
Validation  testing  with  the  AGARD  445.6  flutter  test  case  showed  that  the  aeroelas¬ 
tic  program  performed  well  predicting  flutter  onset  velocity  and  frequency.  Finally, 
initial  testing  of  the  Goland+  LCO  test  case  showed  that  the  program  correctly 
predicted  LCO  behavior. 


3-13 


IV.  Analysis  of  GolandA  LCO 


4-1  Computational  Experiment  Setup 

The  aeroclastic  program  was  used  to  analyze  the  LCO  of  the  Goland+  wing 
described  in  section  3.2.2.  Four  inviscid  grids  and  one  viscous  grid  were  built  to 
look  at  grid  convergence  (Appendix  C)  for  both  steady  solutions  at  angles-of-attack 
from  0  to  20  degrees,  and  for  unsteady  LCO  solutions.  Based  on  the  results  of  this 
study,  it  was  determined  that  the  computational  expense  was  prohibitive  for  viscous 
LCO  simulations.  All  Goland+  LCO  simulations  were  therefore  computed  using  the 
inviscid  solver.  Based  on  the  results  of  the  grid  convergence  study,  the  clean  Goland+ 
results  presented  in  section  4.2  were  obtained  from  an  inviscid  grid  which  consisted 
of  194,780  tetrahedral  cells  with  9,178  cell  faces  on  the  wing. 

At  600  ft/sec  (182.88  m/sec)  and  Mach  0.92,  the  dynamic  pressure  was  be¬ 
yond  the  flutter  point  and  was  shown  by  Snyder,  et  al.  [51]  to  result  in  LCO. 
For  the  Goland+,  these  flow  conditions  fell  within  the  transonic  flutter  dip  region. 
To  match  the  conditions  of  Snyder,  et  al.  [51],  the  far-held  density,  temperature, 
and  pressure  were  set  to  0.0023771  slugs/ ft3  (1.225  kg/m3),  518.67 °R  (288.15  K), 
and  722.1813  lb/ ft2  (34578.04  Pa)  respectively,  while  R  was  set  to  585.7438 
(97.951  f/'f/f ) .  These  how  parameters  were  used  for  all  the  clean  wing  and  store 
simulations  in  this  study. 

All  of  the  Goland+  simulations  used  a  time  step  of  0.001  seconds  corresponding 
to  approximately  300  time  steps  per  cycle.  To  check  for  time-step  convergence,  a 
time  step  of  0.0005  was  tested.  Cutting  the  time  step  in  half  changed  the  mode  1 
response  by  0.367%,  the  mode  2  response  by  0.027%,  the  Cl  response  by  0.029%, 
and  the  frequency  by  0.160%.  These  changes  were  sufficiently  small  that  the  time 
step  of  0.001  seconds  was  considered  converged. 

To  determine  the  role  of  store  aerodynamics  in  LCO,  grids  for  various  store 
configurations  were  created  and  used  with  the  same  structural  model  used  by  the 
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clean  wing.  Therefore,  in  these  cases,  the  stores  were  modelled  aerodynamically  but 
not  structurally.  The  tip  mass  was  retained  in  the  structural  model  because  it  was 
necessary  for  the  clean  wing  to  obtain  LCO.  The  grids  were  created  by  starting  with 
the  clean  grid  of  68,949  tetrahedral  cells  with  518  cell  faces  on  the  wing,  and  then 
adding  store  shapes. 

Under- wing  and  tip  stores  were  added  to  the  Goland+  wing  and  the  same  initial 
perturbation  and  boundary  conditions  were  used  as  for  the  clean  wing.  Under- wing 
stores  can  be  added  to  the  Goland+  wing  at  any  rib  station.  These  stations  are  evenly 
spaced  2  ft  (0.6096  m)  apart.  Aerodynamic  forces  and  moments  were  calculated  for 
the  store  at  the  point  the  pylon  joined  the  wing  on  the  elastic  axis.  These  forces 
and  moments  were  transferred  into  the  structural  model  at  this  point.  Under-wing 
stores  were  added  at  20%  (4  ft  (1.2192  m)),  50%  (10  ft  (3.048  m)),  and  80%  (16  ft 
(4.8768  m))  half-span.  The  under- wing  configurations  consisted  of  the  Goland+ 
grid  with  a  10  ft  (3.048  m)  long,  12  in  (0.3048  m)  diameter  cylindrical  store  with 
an  elliptic  nose  cone  centered  below  the  wing.  The  top  of  the  store  was  located 
12  in  (0.3048  m)  below  the  bottom  of  the  wing  and  was  attached  to  the  wing  via  a 
biconvex-shaped  pylon  that  was  3  ft  (0.9144  m)  long  centered  chordwise  under  the 
wing.  The  grid  with  a  tip  store  consisted  of  the  Goland+  grid  with  a  10  ft  (3.048  m) 
long,  5  in  (0.127  m)  diameter  cylindrical  store  with  an  elliptic  nose  cone  centered  on 
the  wing  tip.  These  four  configurations  are  shown  in  Figure  4.1. 

Additional  store  configurations  were  also  analyzed  in  order  to  determine  how 
store  modifications  change  the  store  forces  or  change  the  interference  on  the  bottom 
of  the  wing.  These  additional  configurations  included  the  addition  of  fins  to  the 
stores,  changes  in  pylon  height,  fore/aft  movement  of  the  stores,  changes  in  store 
diameter,  and  multiple  stores.  The  effect  fins  on  stores  had  on  LCO  was  examined  by 
adding  fins  to  the  tail  of  the  stores  for  the  four  store  configurations.  The  pylon  height 
effect  was  studied  by  changing  the  height  from  12  in  (0.3048  m)  to  6  in  (0.1524  m). 
The  effect  of  stream-wise  positioning  of  the  store  was  examined  by  modifying  the 
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(a)  Goland+  20%  Half-span  Under¬ 
wing  Store 


(b)  Goland+  50%  Half-span  Under¬ 
wing  Store 


Wing  Tip 


Wing  Tip 


Wing  Root 


Wing  Root 


(c)  Goland+  80%  Half-span  Under-  (d)  Goland+  Tip  Store 

wing  Store 

Figure  4.1  Goland+  Wing  with  Aerodynamic  Stores  Surface  Grids 
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Figure  4.2  Goland+  with  Multiple  Stores  Surface  Grids 

position  of  the  under-wing  store  grids  by  shifting  the  stores  2  ft  (0.6096  m)  fore 
and  aft,  while  the  12  in  (0.3048  m)  pylon  remained  centered  beneath  the  wing.  The 
effect  store  diameter  had  on  LCO  was  examined  by  doubling  the  under-wing  store 
diameter  to  24  in  (0.6096  m)  and  the  tip  store  diameter  to  10  in  (0.254  m).  Finally, 
multiple,  small- diameter  stores  were  added  to  the  Goland+  wing  to  create  further 
interference  on  the  bottom  of  the  wing.  The  tip  store  and  the  under-wing  stores  at 
20%  and  50%  half-span,  were  added  to  the  Goland+  wing  (Figure  4.2a).  This  was 
repeated  for  the  large  diameter  stores  (Figure  4.2b). 

Inertial  forces  were  added  to  the  stores  in  order  to  determine  the  relative  im¬ 
portance  of  store  aerodynamics  versus  store  mass.  The  inertial  forces  were  calculated 
based  on  a  mass  of  25  slugs  (364.8476  kg)  located  a  fixed  distance  from  the  elastic 
axis.  The  configurations  with  a  12  in  (0.3048  m)  diameter  under-wing  store  centered 
on  a  12  in  (0.3048  m)  pylon  and  the  5  in  (0.127  m)  diameter  tip  store  were  run  with 
a  store  mass  located  -2  ft  (-0.6096  m),  -1  ft  (-0.3048  m),  0  ft  (0.0  m),  1  ft  (0.3048  m), 
and  2  ft  (0.6096  m)  in  the  x-direction  (stream-wise  direction)  from  the  elastic  axis. 
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4-2  Clean  GolandC  Wing  Results 

The  LCO  for  the  Goland+  clean-wing  configuration  consisted  of  a  nearly  in- 
phase  coupling  of  mode  1  and  mode  2  which  produced  a  single-degree-of-freedom, 
torsional  motion,  with  an  axis  of  rotation  running  from  the  elastic  axis  at  the  wing 
root  to  a  point  slightly  forward  of  the  leading  edge  at  the  wing  tip.  The  LCO  had  an 
amplitude  of  ±1.065  for  Cl  (Figure  4.3),  ±0.47  for  Cm,  and  ±20.6  degrees  for  the  tip 
angle-of-attack.  There  was  significant  hysteresis  in  Cm  as  shown  in  Figure  4.5.  This 
was  expected  because  LCO  is  a  non-linear  phenomenon.  If  it  were  linear  with  no 
hysteresis,  Cm  versus  angle-of-attack  would  form  a  straight  line.  The  two  dominant 
modes  had  an  amplitude  of  ±8.8  for  mode  1  and  ±7.5  for  mode  2  (Figure  4.3). 
The  amplitude  of  mode  3  was  ±0.66  which  was  only  7.5%  of  the  contribution  of 
mode  1.  The  remaining  modes  contributed  even  less.  A  power-spectral-density  plot 
of  the  motion  showed  that  the  power  was  primarily  at  3.3  Hz  and  9.8  Hz  with  some 
contribution  from  6.5  Hz  and  10.7  Hz  (Figure  4.4).  The  dominant  modes  in  the 
LCO  response  were  mode  1  followed  closely  by  mode  2.  However,  the  dominant 
forces  were  mode  1  followed  by  mode  4.  The  remaining  force  contribution  was  in 
mode  3,  mode  5,  then  finally  mode  2  and  mode  6  (Figure  4.6). 

The  only  nonlinearity  present  was  due  to  the  existence  and  movement  of  shocks. 
There  was  no  separation  (inviscid  flow)  and  no  structural  nonlinearity.  The  shock 
was  initially  located  at  approximately  80%  chord  where  it  ran  from  the  wing  root, 
to  the  wing  tip,  with  a  slight  curve  forward  before  disappearing  at  the  wing  tip 
(Figure  4.7).  As  the  wing  began  to  deform,  the  shock  started  to  slowly  move  fore 
and  aft  sinusoidally  (type-A  shock  motion  [56]).  As  the  amplitude  of  the  wing  motion 
increased,  the  magnitude  of  the  shock  strength  changed  periodically  with  the  motion 
until  the  shock  disappeared  during  part  of  its  aft  motion  (type-B  shock  motion  [56]). 
Transition  from  type-A  to  type-B  was  observed  at  approximately  2.9  seconds,  ranging 
from  the  wing  root  to  about  50%  span.  The  shock  motion  for  the  outboard  50% 
span  remained  type-A.  The  inboard  type-B  shock  motion  limited  the  energy  flow 
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Goland+  Clean  Wing 
Mach  0.92,  600  ft/sec 


Figure  4.3  Clean  Wing  Cl  and  Modal  Amplitudes 


PSD  of  CL,  0.0  to  10.0  sec 


Figure  4.4  Clean  Wing  PSD 
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Figure  4.5  Clean  Wing  Cm  vs  Angle-of-Attack 


Figure  4.6  Clean  Wing  Modal  Force  Amplitudes 
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Figure  4.7  Steady  State  Density  Contours 

and  changed  the  mode  1  amplitude  growth  from  an  exponential  growth  to  a  linear 
growth,  as  seen  in  Figure  4.3.  This  type-B  motion  continued  to  gain  strength  until 
it  consisted  of  an  appearing  and  disappearing  shock  at  the  trailing  edge  with  no 
observable  fore/aft  motion.  However,  this  type-B  shock  motion  was  only  on  the 
inboard  50%  span.  Most  of  the  energy  flowing  into  the  structure  was  occurring  on  the 
outboard  50%  span.  At  approximately  4.6  seconds,  a  periodic  lambda  shock  began 
to  appear  and  disappear  near  the  wing  tip  in  conjunction  with  a  strong  normal  shock 
moving  fore  and  aft  at  the  trailing  edge.  This  lambda  shock  and  normal  shock  further 
limited  the  energy  transferred  from  the  fluid  to  the  structure.  At  6.4  seconds,  the 
wing  was  in  LCO.  The  combination  of  type-B  shock  motion,  lambda  shock  motion, 
and  trailing-edge  normal-shock  motion  were  responsible  for  quenching  the  energy 
flow  from  the  fluid  to  the  structure  resulting  in  LCO  for  this  wing. 

The  coupling  of  mode  1  and  mode  2  into  a  single-degree-of-freedom,  torsional 
motion  about  a  point  ahead  of  the  leading  edge,  leads  to  a  wing  twisting  motion  such 
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that  the  aerodynamic  forces  on  the  wing  are  constantly  pushing  the  wing  toward  the 
static  aeroelastic  position.  When  the  wing  bends  down  and  twists  up,  the  aerody¬ 
namic  forces  try  to  restore  the  wing  to  its  static  aeroelastic  position.  However,  the 
restoring  forces  are  large  enough  that  the  amount  of  energy  transferred  into  the  wing 
causes  the  wing  to  overshoot  its  neutral  point,  bending  and  twisting  even  further  in 
the  opposite  direction,  leading  to  a  growing  oscillation.  When  the  oscillations  grow 
to  the  point  that  the  trailing-edge  shock  and  lambda  shock  appear  near  the  wing 
tip,  the  shocks  limit  the  force  on  the  wing,  thereby  limiting  the  overshoot.  When  the 
shock  strength  reaches  the  point  that  it  balances  the  inertial  forces,  a  stable  LCO  is 
achieved.  Figure  4.8  illustrates  the  pressure  distribution  on  the  wing.  If  the  normal 
shock  and  lambda  shock  were  not  present,  the  coefficient  of  pressure  (Cp)  on  the  top 
of  the  wing  would  be  lower.  The  net  force  on  the  wing  would  therefore  be  higher, 
leading  to  the  wing  bending  and  twisting  further  in  each  subsequent  cycle. 

During  a  LCO  cycle  (Figure  4.9),  the  wing  starts  bending  down  and  twisting 
up.  A  normal  shock  appears  near  the  trailing  edge  on  the  outboard  portion  of  the 
wing  (Figure  4.9a).  This  normal  shock  increases  the  pressure  on  top  of  the  wing  and 
begins  to  limit  the  flow  of  energy  from  the  fluid  into  the  wing.  As  the  wing  continues 
to  bend  down  and  twist  up,  a  lambda  shock  appears  near  the  tip  (Figure  4.9b).  When 
the  lambda  shock  first  appears,  it  runs  from  the  leading  edge  of  the  wing  near  the 
tip  to  where  the  trailing-edge  normal  shock  meets  the  trailing  edge  near  50%  span. 
As  the  wing  continues  to  bend  down  and  twist  up,  the  lambda  shock  strengthens 
and  the  shock  angle  increases.  At  the  extreme  of  the  wing  motion,  the  lambda  shock 
runs  from  the  leading  edge  of  the  wing  near  the  tip,  to  the  center  of  the  normal 
shock  at  the  trailing-edge,  at  approximately  75%  span  (Figure  4.9c).  As  the  wing 
then  begins  to  straighten  and  untwist,  the  lambda  shock  reverses  direction  and  the 
shock  angle  begins  to  decrease.  The  normal  shock  at  the  trailing  edge  also  begins 
to  weaken  and  move  aft  (Figure  4.9d).  At  the  end  of  the  the  first  half  of  the  LCO 
cycle,  the  shocks  are  completely  gone  and  the  forces  on  the  top  and  bottom  of  the 
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Figure  4.8  Clean  Wing  at  6.815  Seconds  (Pressure  Distribution) 


wing  are  nearly  balanced  (Figure  4.9e).  This  is  repeated  for  the  second  half  of  the 
LCO  cycle  when  the  wing  starts  to  bend  up  and  twist  down. 

Figure  4.10  shows  the  steady  pressure  contours  for  a  rigid  wing  at  20  degrees 
angle-of-attack.  The  forces  on  the  rigid-wing  airfoil  section  are  much  higher  than 
those  shown  in  Figure  4.8e,  even  though  the  angle-of-attack  at  the  tip  is  approxi¬ 
mately  the  same.  This  difference  in  force  is  because  of  the  shocks,  which  provide 
the  quenching  mechanisms  responsible  for  LCO  in  this  wing.  The  normal  shock  and 
lambda  shock  do  not  appear  if  the  wing  is  given  a  constant  angle-of-attack.  However, 
if  the  wing  is  deformed  into  the  bent  and  twisted  shape  that  occurs  during  LCO  and 
a  steady-state  flow  solution  is  computed,  the  lambda  and  trailing-edge  shocks  are 
present.  These  shocks  are  present  because  of  the  shape  of  the  wing  and  not  because 
of  the  dynamic  motion  of  the  wing.  The  dynamic  motion  does,  however,  play  a  role 
in  the  shock  strength  and  movement.  The  shocks  are  weaker  and  slightly  aft  when 
computed  for  steady-state  flow  versus  those  computed  for  dynamic  motion.  The 
coupling  of  the  primary  bending  and  twisting  modes  is  what  leads  to  this  deformed 
shape,  which  in  turn  leads  to  shock  motion  that  provides  the  energy  quenching  nec¬ 
essary  for  LCO.  The  aerodynamic  forces  dictate  how  fast  the  oscillations  grow  and 
the  final  amplitudes,  not  whether  LCO  exists.  This  was  further  demonstrated  by 
changing  the  airfoil  on  the  Goland+  wing  in  order  to  change  the  aerodynamic  forces 
(Appendix  D). 

If  a  wing  is  undergoing  LCO,  a  phase  plot  will  trace  an  orbit  along  a  closed  path 
and  there  will  be  zero  net  work  over  a  cycle  [5].  For  the  clean- wing  configuration, 
plunge  rate  versus  plunge  was  plotted  for  the  point  where  the  elastic  axis  intersects 
the  wing  tip.  This  plunge  phase  plot,  shown  in  Figure  4.11a,  shows  the  growth  of  the 
amplitudes  culminating  in  a  closed  path  indicating  LCO.  The  pitch  rate  of  the  wing 
tip  was  also  plotted  versus  pitch.  This  pitch  phase  plot,  shown  in  Figure  4.11b,  also 
shows  the  growth  of  the  amplitudes  culminating  in  LCO.  Because  the  LCO  motion 
appears  to  be  a  single-degree-of-freedom,  torsional  motion,  the  work  being  done  on 
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Figure  4.9  Clean  Wing  (Density  Contours) 
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20.0  Degrees  Angle  of  Attack 


Chord 

Figure  4.10  Rigid  Wing  Airfoil  Cp  plot  for  20  degrees  Angle-of- Attack 

the  wing  by  the  fluid  can  be  compared  with  the  wing-tip  angle-of-attack.  As  can  be 
seen  in  Figure  4.12,  whenever  the  wing  is  twisting  away  from  zero  degrees,  work  is 
being  done  on  the  fluid,  and  whenever  the  wing  is  twisting  toward  zero  degrees,  work 
is  being  done  on  the  structure.  This  results  in  two  work  cycles  per  LCO  cycle.  There 
is  a  5-degree  phase  lag  between  the  work  and  angle-of-attack.  This  is  because  the 
bending  mode  (mode  1)  lags  the  torsion  mode  (mode  2)  by  5  degrees  (Figure  4.3). 
Integrating  the  work  over  an  LCO  cycle  should  result  in  zero  net  work.  However, 
for  this  case,  the  work  did  not  sum  to  zero.  There  was  an  error  of  0.05%  of  the  total 
work  done.  This  was  because  the  thin-plate  spline  did  not  conserve  energy  across 
the  grid  boundaries. 

4-3  Goland+  Wing  with  Under-Wing  Stores  Results 

One  of  the  goals  of  this  study  was  to  determine  the  role  of  store  aerodynamics 
in  LCO.  It  was  found  that  aerodynamic  store  shapes  affect  LCO  in  two  offsetting 
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Goland+  Clean  Wing 
Mach  0.92,  600  ft/sec 


(b)  Pitch  Phase  Plot  (Wing  Tip) 

Figure  4.11  Clean  Wing  Phase  Plots 
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Time 


(a)  Work  ( N  ■  m)  on  Structure  for  10  seconds 


Time 


(b)  Work  ( N  ■  m)  on  Structure  for  One  LCO  Cycle  (Expanded 
Region) 

Figure  4.12  Work  on  Goland+  Structure  (Clean  Wing) 

4-15 


Table  4.1  Effect  of  Under-Wing  Stores  on  Peak  Values 


20%  half-span 

50%  half-span 

80%  half-span 

CL 

1.10% 

4.08% 

1.28% 

Tip  AOA 

1.12% 

4.96% 

2.56% 

Mode  1  Force 

1.48% 

5.22% 

0.28% 

Mode  2  Force 

13.69% 

25.15% 

14.93% 

Mode  3  Force 

2.44% 

1.01% 

1.73% 

Mode  4  Force 

0.50% 

1.96% 

0.57% 

ways:  by  interfering  with  the  flow  field  on  the  wing  surface,  and  by  transferring 
additional  store  forces  into  the  structure.  The  under-wing  stores  interfere  with  the 
airflow  on  the  lower  surface  of  the  wing  which  decreases  LCO  amplitudes,  whereas, 
store  forces  transferred  into  the  wing  structure  directly  increase  LCO  amplitudes. 

The  further  outboard  the  store  was  placed,  the  greater  the  angle-of-attack 
experienced  by  the  store.  This  resulted  in  higher  store  loads  and  increased  LCO  am¬ 
plitude.  An  exception  to  this  was  found  at  the  50%  half-span  location  (Figure  4.13). 
Small  stores  with  the  12  in  (0.3048  m)  pylon  were  placed  at  20%,  50%,  and  80% 
half-span  locations  as  discussed  in  section  4.1.  Adding  stores  increased  the  mode  2 
response  in  all  cases.  The  mode  1  response  decreased  for  stores  at  20%  and  80% 
half-span  and  increased  for  stores  at  50%.  It  was  speculated  that  this  increase  was 
because  of  the  mode  3  mode-shape  (Figure  3.7)  which  had  large  displacements  in  the 
middle  of  the  wing.  Further  weight  was  added  to  this  speculation  by  observing  that 
the  stores  primarily  contributed  to  mode  3  forces  implying  that  the  store  located 
where  mode  3  provided  the  greatest  deformation  had  the  greatest  effect.  Adding 
the  under-wing  stores  increased  the  magnitude  of  the  total  modal  forces,  the  magni¬ 
tude  of  Cl,  and  the  magnitude  of  the  tip  angle-of-attack  (Table  4.1).  These  changes 
resulted  in  LCO  being  reached  more  quickly  than  with  a  clean  wing. 

For  the  20%  and  80%  half-span  configurations,  the  stores  only  contributed 
0.2%  of  the  peak  force  in  mode  3.  The  50%  store  contributed  0.7%  of  the  peak 
force  in  mode  3.  To  increase  these  store  forces,  fins  were  added  to  the  tail  of  all 
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the  store  configurations.  The  fins  were  aft  of  the  wing  so  they  did  not  affect  the 
flow  field  on  the  wing.  For  these  small  stores,  fins  did  not  make  a  difference  in  the 
LCO  behavior  (Figure  4.14).  They  did  increase  the  store  forces,  but  the  store  forces 
were  so  small  compared  to  the  wing  forces  that  the  effects  on  LCO  were  negligible. 
The  pylon  height  was  changed  from  12  in  (0.3048  m)  to  6  in  (0.1524  m)  in  order  to 
change  the  store  moments.  Shortening  the  pylon  only  made  a  small  difference  in  the 
LCO  behavior  (Figure  4.15).  The  store  forces  and  the  wing  response  were  similar, 
regardless  of  the  pylon  height.  Shortening  the  pylon  led  to  lower  mode  2  amplitudes, 
or  less  twist,  but  otherwise,  the  behavior  was  the  same  regardless  of  pylon  height. 
The  fins  and  pylon  height  weakly  affected  the  twisting  moment,  whereas,  the  bending 
force  was  the  primary  contributor  to  LCO. 

The  store  forces  that  were  transferred  into  the  structural  model  were  magnified 
in  order  to  amplify  how  the  store  forces  affected  the  LCO.  Solutions  for  the  wing  with 
a  12  in  (0.3048  m)  diameter  under-wing  store  centered  on  a  12  in  (0.3048  m)  pylon 
were  computed,  but  with  the  store  forces  that  were  calculated  from  the  aerodynamic 
loads  being  magnified  by  ±5  before  being  transferred  into  the  structural  model. 
Altering  the  store  forces  in  this  manner  affected  the  LCO  (figures  4.16  and  4.17).  The 
store  forces  increased  the  total  bending  force  on  the  wing,  which  increased  the  LCO 
amplitude.  The  bending  force  primarily  affected  LCO  by  directly  contributing  to 
mode  1.  Magnifying  the  bending  force  increased  the  amplitude  of  the  LCO  because 
it  increased  the  restoring  force,  resulting  in  greater  overshoot  beyond  the  neutral 
position.  Magnifying  and  reversing  the  phase  of  the  loads  (multiplying  forces  by  —5) 
decreased  the  restoring  force  resulting  in  less  overshoot  and  a  decreased  amplitude 
LCO.  These  results  reinforced  the  hypothesis  that  the  increased  LCO  amplitudes 
observed  with  small  under-wing  stores  were  caused  by  the  store  loads  which  were 
greater  than  the  quenching  effect  caused  by  interference  of  the  flow  on  the  bottom 
of  the  wing. 
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(b)  Pitch  Phase  Plot  (Wing  Tip) 

Figure  4.14  Phase  Plots  for  Wing  with  Finned-Stores 
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(a)  Plunge  Phase  Plot  (Wing  Tip) 
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(b)  Pitch  Phase  Plot  (Wing  Tip) 


Figure  4.15  Phase  Plots  for  Pylon  Height  Chang 
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(a)  Plunge  Phase  Plot  (Wing  Tip) 


(b)  Pitch  Phase  Plot  (Wing  Tip) 


Figure  4.16  Phase  Plots  when  Store  Forces  Magnified  by  -5 
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The  second  mechanism  by  which  under-wing  stores  affected  LCO  was  by  inter¬ 
fering  with  the  airflow  on  the  bottom  of  the  wing.  The  80%  half-span  configuration 
with  a  small  store  was  computed  with  the  store  forces  turned  off  in  order  to  isolate 
this  interference  effect.  When  the  store  forces  were  ignored,  Cl  decreased  0.97%, 
Cm  decreased  1.26%,  the  tip  angle-of-attack  decreased  1.59%,  mode  1  amplitudes 
decreased  1.67%,  and  mode  2  amplitudes  decreased  1.82%.  These  results  demon¬ 
strated  that  interference  from  a  small  store  slightly  quenched  the  LCO. 

The  under-wing  store  grids  with  the  12  in  (0.3048  m)  pylon  and  fins  were 
modified  to  see  how  stream- wise  store  position  affects  the  interference  on  the  bottom 
of  the  wing.  The  stores  were  shifted  2  ft  (0.6096  m)  fore  and  aft.  The  pylon  remained 
centered  beneath  the  wing.  Moving  the  store  aft  produced  LCO  (Figure  4.18).  When 
the  50%  half-span  store  and  the  80%  half-span  store  were  moved  forward,  the  flow 
separated  on  the  store  causing  the  inviscid  fluid  solver  to  fail,  so  no  conclusions  could 
be  drawn.  The  20%  half-span  store  resulted  in  LCO  when  moved  forward.  This  LCO 
magnitude  was  similar  to  the  original  centered  store  position.  When  the  stores  were 
moved  aft,  they  only  slightly  affected  the  LCO.  Shifting  under-wing  stores  fore  or 
aft  only  slightly  changed  the  LCO  because  the  change  in  the  forces  was  very  small 
and  the  interference  of  the  flow  on  the  wing  because  of  the  store  remained  small. 

To  further  amplify  the  interference  effect,  the  stores  without  fins  were  doubled 
in  diameter  to  24  in  (0.6096  m).  The  centered  12  in  (0.3048  m)  pylon  was  retained. 
As  can  be  seen  in  Figure  4.19,  the  large  diameter  store  had  a  large  affect  on  the 
LCO.  Whereas  adding  a  small  under-wing  store  produced  an  LCO  of  slightly  greater 
amplitude,  adding  a  large  under-wing  store  resulted  in  an  LCO  of  smaller  amplitude. 
The  further  outboard  the  large  store,  the  smaller  the  amplitude  of  the  LCO.  Again, 
the  total  modal  forces  did  not  change  greatly  for  large  stores.  The  aerodynamic 
forces  were  slightly  higher  with  a  large  store  than  with  a  small  store  but  were  still 
less  than  1%  of  the  total  forces  for  all  modes.  This  indicated  that  the  large  store  was 
affecting  the  flow  on  the  bottom  of  the  wing  and  it  was  influencing  the  dynamic  shock 
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(a)  Plunge  Phase  Plot  (Wing  Tip) 


Pitch 


(b)  Pitch  Phase  Plot  (Wing  Tip) 

Phase  Plots  for  Wing  with  Stores  Shifted  Stream-Wise 


motion  which  was  limiting  the  energy  flow  into  the  structure.  The  interference  effect 
was  not  great  enough  to  damp  out  the  oscillations,  but  it  was  enough  to  decrease 
the  LCO  amplitude. 

The  large  store  at  80%  half-span  illustrates  how  under-wing  stores  affect  the 
flow  field  on  the  bottom  of  the  wing,  and  thereby,  affect  the  LCO.  When  the  wing 
was  bent  up  and  twisted  down,  the  flow  on  the  top  of  the  wing  was  unchanged  from 
the  clean  wing  condition.  However,  on  the  bottom  of  the  wing  (Figure  4.20a)  there 
were  significant  changes  when  compared  to  the  clean  wing  at  the  same  point  in  an 
LCO  cycle  (Figure  4.20b).  The  lambda  shock  was  weak  if  detectable  at  all,  and  the 
normal  shock  at  the  trailing  edge  was  weaker  than  that  experienced  by  the  clean 
wing.  The  store  wake,  or  the  flow  from  the  store  that  impinged  on  the  flow  on  the 
bottom  of  the  wing,  increased  the  pressure  on  the  front  half  of  the  wing  in  the  vicinity 
of  the  store.  The  normal  shock  then  increased  the  pressure  again  (Figure  4.21).  The 
increase  in  pressure  decreased  the  downward  restoring  force  on  the  wing,  which  had 
the  same  effect  as  the  lambda  shock  on  the  clean  wing,  limiting  the  energy  flow 
into  the  wing.  When  the  wing  was  bent  down  and  twisted  up,  the  flow  on  the  top 
of  the  wing  was  unchanged  from  the  clean  wing  condition.  The  lambda  shock  and 
trailing-edge  shock  were  present,  though  weaker  than  those  present  in  the  clean  wing 
because  of  the  lower  pitch  and  plunge  amplitudes  (Figure  4.22).  These  shocks  were 
still  providing  quenching  leading  to  LCO.  If  the  store  was  large  enough  to  have  a 
large  wake  raise  the  pressure  on  the  wing,  it  provided  quenching.  A  very  small  store 
with  a  minimal  wake  had  minimal  impact  on  the  LCO  caused  by  interference  from 
the  store. 

Multiple,  small-diameter  stores  were  added  to  the  Goland+  wing  to  create 
further  interference.  A  tip  store  and  the  under-wing  stores  at  20%  and  50%  half¬ 
span  were  added  to  the  Goland+  wing.  The  centered  12  in  (0.3048  m)  pylon  was 
retained.  The  small-diameter  stores  resulted  in  an  LCO  with  an  amplitude  slightly 
greater  than  the  clean  wing.  For  the  small  stores,  the  store  forces  effect  was  stronger 
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Large  80%  Store  at  10.495  Seconds  (Density  Contours) 
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(e)  5.5m  span 


Figure  4.21 
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Figure  4.22 
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than  the  interference  damping  effect  (Figure  4.23).  Multiple,  large-diameter  stores 
were  then  added  to  the  Goland+  wing  to  create  further  interference.  When  the  large- 
store  configuration  was  run,  the  oscillations  damped  out.  The  multiple,  large  stores 
interfered  with  the  flow  on  the  bottom  of  the  wing  to  the  extent  that  it  prevented 
the  energy  flow  into  the  structure,  resulting  in  damping. 

One  effect  that  under-wing  stores  had  on  LCO  was  shown  by  the  mode  1 
response  oscillating  about  a  negative  value  instead  of  about  zero,  as  with  a  clean 
wing  or  with  a  tip  store.  This  showed  that  the  presence  of  the  under-wing  stores 
shifted  the  static  aeroelastic  solution.  Shifting  the  static  aeroelastic  solution  did 
not  change  the  LCO  response,  only  the  neutral  configuration  about  which  the  wing 
oscillated.  This  was  further  demonstrated  by  computing  the  solution  for  the  clean 
Goland+  wing  at  a  positive  angle-of- attack  (Appendix  E). 

To  examine  the  effect  of  store  mass  position  on  LCO  as  opposed  to  the  store 
aerodynamics,  inertial  forces  were  added  to  the  store  aerodynamic  forces.  The  12  in 
(0.3048  m)  diameter  under-wing  stores  centered  on  12  in  (0.3048  m)  pylons  were 
run  with  a  store  mass  located  -2  ft  (-0.6096  m),  -1  ft  (-0.3048  nr),  0  ft  (0.0  nr),  1  ft 
(0.3048  nr),  and  2  ft  (0.6096  nr)  in  the  x-direction  (stream-wise  direction)  from  the 
elastic  axis.  As  can  be  seen  in  Figure  4.24,  when  the  mass  was  added  on  the  elastic 
axis  or  forward  of  the  elastic  axis,  the  LCO  amplitude  only  slightly  increased  for 
the  under-wing  store  located  at  20%  half-span.  However,  if  the  mass  was  added  1  ft 
(0.3048  nr)  aft  of  the  elastic  axis,  a  different  LCO  mode  was  encountered  (figures  4.24 
and  4.25).  The  amplitudes  grew  in  the  same  manner  as  the  case  without  the  mass, 
but  at  approximately  4  seconds,  mode  4  started  growing  and  eventually  stabilized. 
At  approximately  14  seconds,  a  new  LCO  state  was  reached.  When  the  mass  was 
added  2  ft  (0.6096  nr)  aft  of  the  elastic  axis,  all  of  the  modes  grew,  leading  to 
divergent  flutter. 

When  the  mass  was  added  to  the  50%  half-span,  under-wing  store,  divergent 
flutter  was  obtained  for  all  stream-wise  locations.  It  was  hypothesized  that  this  was 
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(a)  Plunge  Phase  Plot  (Wing  Tip) 


- B -  Clean  Wing 

Goland+  Wing  - A -  Store20  Base 

Mach  0.92,  600  ft/sec  - V —  Store20  2’  Mfore 


Pitch 


(b)  Pitch  Phase  Plot  (Wing  Tip) 

Figure  4.24  Phase  Plots  for  Store  Point  Mass  at  20%  Half-span 
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Figure  4.25  Amplitudes  for  Point  Mass  1’  Aft  of  EA  at  20%  Half-span 

because  the  stores  primarily  contributed  to  mode  3  forces.  When  the  store  mass  was 
located  where  mode  3  provided  the  greatest  contribution,  it  had  a  large  effect.  When 
the  mass  was  placed  on  the  elastic  axis,  it  came  close  to  stabilizing  into  an  LCO. 
It  initially  appeared  to  enter  the  LCO  condition  observed  without  the  mass,  but  at 
approximately  6  seconds,  mode  4  started  to  diverge  followed  by  mode  3,  eventually 
driving  the  system  to  divergent  flutter  (Figure  4.26).  For  the  other  50%  half-span 
mass  locations,  the  system  immediately  went  into  divergent  flutter. 

Adding  the  mass  to  the  80%  half-span  store  illustrated  the  stabilizing  and 
destabilizing  effect  mass  position  provided,  as  shown  in  Figure  4.27.  When  the  mass 
was  located  on  the  elastic  axis,  a  LCO  response  similar  to  the  no-mass  case  was 
observed.  As  the  mass  was  moved  forward  of  the  elastic  axis,  the  LCO  amplitudes 
decreased  drastically  while  mode  2  increased  in  dominance  and  mode  1  decreased. 
This  result  was  expected  because  adding  a  mass  forward  of  the  elastic  axis  raised  the 
flutter  speed,  which  led  to  a  decreased  amplitude  LCO.  The  modes  continued  to  be 
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Goland+  50%  Store,  12"  Pylon,  Mass  on  EA 
Mach  0.92,  600  ft/sec 


Figure  4.26  Amplitudes  for  Point  Mass  on  EA  at  50%  Half-span 

coupled  into  a  single-degree-of-freedom  response.  When  the  mass  was  moved  aft  of 
the  elastic  axis,  the  LCO  amplitudes  increased  significantly.  Again,  this  matched  the 
expected  result  because  adding  a  mass  aft  of  the  elastic  axis  lowered  the  flutter  speed, 
leading  to  an  increased  amplitude  LCO.  Mode  1  continued  to  grow  in  dominance 
while  mode  2  decreased  in  importance.  The  modes  continued  to  be  coupled  into  a 
single-degree-of-freedom  response,  but  with  mode  4  growing  and  beginning  to  play 
a  role  in  the  response.  However,  by  the  time  the  mass  had  moved  2  ft  (0.6096  m) 
aft,  it  was  no  longer  a  single-degree-of-freedom  response,  and  a  different  LCO  state 
was  reached.  If  the  mass  continued  to  move  aft,  mode  4  would  begin  to  dominate 
and  divergent  flutter  would  be  encountered.  Adding  inertial  forces  did  not  change 
the  shock  motion  nor  change  the  quenching  mechanism.  It  did  change  the  inertia  of 
the  wing  and  determined  how  efficient  the  shock  motion  was  at  limiting  the  energy 
transferred  into  the  structure.  As  the  wing  twisted  more  due  to  a  forward  mass, 
the  shock  motion  caused  the  energy  flow  to  be  limited  before  the  inertia  grew  too 
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high,  decreasing  the  amplitude  of  the  LCO.  Conversely,  as  the  wing  bent  more  for  a 
rearward  mass,  the  inertia  grew  higher  before  the  shock  motion  limited  the  energy 
flow,  thereby,  increasing  the  amplitude  of  the  LCO. 

Jh4  Goland+  Wing  with  Tip- Store  Results 

A  tip  store  affected  the  circulation  of  the  fluid  around  the  wing  tip  and  con¬ 
tributed  to  store  forces,  but  it  did  not  interfere  with  the  airflow  on  the  bottom  of 
the  wing.  Adding  a  tip  store  increased  the  mode  1  and  mode  2  responses.  Like  the 
under-wing  stores,  the  tip  store  primarily  contributed  to  force  in  mode  3  where  the 
tip  store  was  responsible  for  4.8%  of  the  total  mode  3  force.  However,  the  total  force 
in  mode  3  decreased  by  6.63%  from  that  of  the  clean  wing.  Adding  the  tip  store 
increased  Cl  by  5.06%  and  the  tip  anglc-of-attaek  by  1.43%.  Unlike  the  under-wing 
store,  the  mode  1  response  for  the  tip  store  oscillates  about  zero  instead  of  about 
a  negative  value,  the  same  as  the  clean  wing.  These  changes  resulted  in  the  am¬ 
plitudes  of  the  motion  growing  faster  than  the  clean  wing  and  reaching  LCO  more 
quickly.  The  tip  store  increased  the  mode  1  amplitudes  which  increased  the  amount 
of  bending. 

The  tip  store  was  shifted  2  ft  (0.6096  m)  fore  and  aft.  When  the  tip  store  was 
moved  forward,  the  motion  damped  out.  Analysis  showed  that  the  flow  around  the 
tip  of  the  wing  was  very  sensitive  to  the  presence  of  a  store.  Moving  the  tip  store 
forward  appeared  to  prevent  the  energy  from  flowing  into  the  wing  because  mode  1 
and  mode  2  did  not  couple  into  a  single-degree-of-freedom  flutter  mode,  but  remained 
independent.  LInlike  the  LCO  cases,  store  forces  also  directly  correlated  with  modal 
response  with  mode  1  and  mode  2  being  the  dominant  force  contributors.  Moving 
the  tip  store  aft  resulted  in  an  LCO  state  with  a  similar  twist  and  greater  plunge 
compared  to  the  centered  store  LCO.  In  this  case,  mode  1  increased  greatly  when 
the  store  was  shifted  aft.  The  tip  angle-of-attack,  CL ,  Cm,  and  the  modal  forces  all 
increased  slightly.  Mode  3  decreased  slightly  while  modes  2,  4,  5,  and  6  all  increased 
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slightly.  Overall,  moving  the  tip  store  fore  and  aft  showed  that  the  wing-tip  region 
was  very  sensitive  to  flow  changes.  These  flow  changes  manifested  as  changes  in  the 
mode  1  response  leading  to  damping  (move  tip  store  forward)  or  by  increasing  the 
amount  of  plunge  (move  tip  store  aft). 

In  order  to  amplify  the  effect  the  store  forces  were  having  on  the  LCO,  the 
store  forces  were  magnified  by  ±5  before  being  transferred  into  the  structural  model. 
Magnifying  the  loads  in  this  manner  affected  the  LCO  (figures  4.16  and  4.17).  Tip 
stores  have  the  same  effects  as  under-wing  stores:  magnifying  the  bending  force  on 
the  wing,  and  increasing  the  LCO  amplitude.  The  same  effect  was  seen  when  the 
tip  store  diameter  was  increased  to  10  in  (0.254  m).  As  can  be  seen  in  Figure  4.19, 
the  large  diameter  tip  store  produced  a  higher  amplitude  LCO.  The  large  tip  store 
contributed  9%  of  mode  3  forces.  Unlike  the  under- wing  stores,  this  was  a  significant 
fraction  of  the  total  forces  on  the  wing.  The  store  basically  acted  as  an  extension 
of  wing  area,  increasing  the  bending  forces  on  the  wing  and  increasing  the  LCO 
amplitude. 

To  examine  the  effect  of  tip-store  mass  position  on  LCO  as  opposed  to  store 
aerodynamics,  inertial  forces  were  added  to  the  store  aerodynamic  forces.  The  5  in 
(0.127  m)  diameter  tip  store  was  run  with  a  store  mass  located  -2  ft  (-0.6096  m), 
-1  ft  (-0.3048  m),  0  ft  (0.0  m),  1  ft  (0.3048  m),  and  2  ft  (0.6096  m)  in  the  x-direction 
(stream-wise  direction)  from  the  clastic  axis.  When  the  mass  was  added  at  the  wing 
tip  on  the  elastic  axis,  mode  6  began  to  grow  at  2.5  seconds  and  led  to  divergent 
flutter.  When  the  mass  was  added  aft  of  the  clastic  axis,  all  of  the  modes,  except 
for  mode  2,  began  to  grow  which  led  to  divergent  flutter.  When  the  mass  was  added 
forward  of  the  clastic  axis  (Figure  4.28),  it  appeared  that  a  coupled  LCO  response 
was  obtained  with  mode  2  being  increasingly  dominant  the  further  forward  the  mass 
was  moved.  However,  at  approximately  10  seconds,  mode  4  rapidly  began  to  grow 
leading  to  divergent  flutter  (Figure  4.29).  Adding  inertial  forces  forward  of  the 
clastic  axis  changed  the  coupling  of  the  modes  and  did  not  lead  to  a  lower  amplitude 
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LCO.  As  seen  with  the  under-wing  stores,  adding  the  inertial  forces  did  not  change 
the  shock  motion  nor  change  the  damping  mechanism  but  instead  changed  how  the 
modes  couple.  The  coupling  of  the  primary  bending  and  twisting  modes  is  what  led 
to  shock  motion  that  provided  the  energy  quenching  necessary  for  LCO. 
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Figure  4.29  Amplitudes  for  Point  Mass  1’  Fore  of  EA  at  Tip 
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V.  Conclusions 


5.1  Conclusions 

FLUENT  6.1  was  successfully  integrated  with  a  modal  structural  model  to 
form  an  aeroelastic  program  based  on  a  fully  unstructured  grid  formulation  capable 
of  simulating  flutter  and  LCO.  Using  this  program,  the  aerodynamic  nonlinearity  re¬ 
sponsible  for  LCO  of  the  Goland+  clean-wing  (with  tip  mass)  at  Mach  0.92,  600  ft/sec 
(182.88  m/sec),  was  analyzed.  The  LCO  consisted  of  a  nearly  in-phase  coupling  of 
mode  1  and  mode  2,  which  produced  a  singlc-degree-of-freedom  torsional  motion 
with  an  axis  of  rotation  running  from  the  elastic  axis  at  the  root  of  the  wing,  to 
a  point  slightly  forward  of  the  leading  edge  at  the  tip  of  the  wing.  For  LCO  to 
exist,  a  nonlinearity  must  have  been  present  to  provide  the  quenching  mechanism 
that  limited  the  amplitude  of  the  oscillations.  For  the  Goland+  wing,  the  only  non¬ 
linearity  present  was  the  appearance/disappearance  of  shocks.  Flow  separation  and 
structural  nonlinearity  were  not  modelled.  Periodically  appearing  and  disappearing 
lambda  and  trailing-edge  shocks  decreased  the  restoring  forces  on  the  Goland+  wing, 
resulting  in  a  balancing  of  the  inertial  forces  and  a  stable  LCO. 

For  the  Goland+  wing,  as  the  oscillations  grew,  there  was  a  transition  from 
Tijdeman  type-A  shock  motion  [56],  to  Tijdeman  type-B  shock  motion  [56]  in  the 
inboard  half  of  the  wing.  The  shock  motion  in  the  outboard  half  of  the  wing  span 
remained  type-A.  The  type-B  shock  motion  limited  the  flow  of  energy  from  the  fluid 
into  the  structure,  and  changed  the  amplitude  growth  from  an  exponential  growth  to 
a  linear  growth.  The  type-B  motion  continued  to  gain  strength  until  it  consisted  of 
an  appearing  and  disappearing  shock  at  the  trailing  edge,  with  no  observable  fore/aft 
motion. 

Two  additional  quenching  mechanisms  later  appeared  in  the  outboard  half¬ 
span,  which  further  limited  the  amplitude  growth  and  led  to  LCO.  The  first  was 
a  strong  trailing-edge  shock  that  moved  fore  and  aft,  and  the  second  was  a  strong 
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lambda  shock  that  ran  from  the  leading  edge  of  the  wing  near  the  tip,  to  the  normal 
shock  on  the  trailing  edge.  These  shocks  decreased  the  restoring  force  on  the  wing 
limiting  the  overshoot,  thereby  balancing  the  inertial  forces  producing  a  stable  LCO. 
The  combination  of  type-B  shock  motion,  lambda  shock  motion,  and  the  motion  of 
the  normal  shock  at  the  trailing  edge  was  responsible  for  limiting  the  energy  flow 
from  the  fluid  to  the  structure  resulting  in  LCO. 

Aerodynamic  store  shapes  were  added  to  the  Goland+  wing  at  different  span- 
wise  locations  to  determine  how  store  aerodynamics  affect  LCO.  It  was  found  that 
aerodynamic  store  shapes  affect  LCO  in  two  ways:  by  interfering  with  the  flow  held 
on  the  wing  surface,  and  by  adding  loads  into  the  structure.  Under-wing  stores  inter¬ 
fered  with  the  airflow  on  the  lower  surface  of  the  wing,  diminishing  LCO  amplitudes. 
On  the  other  hand,  the  transfer  of  store  loads  into  the  wing  structure  increased  LCO 
amplitudes  for  both  tip  stores  and  under-wing  stores.  Under-wing  stores  also  led  to  a 
negative  offset  of  the  amplitudes  of  the  primary  bending  mode  by  shifting  the  static 
aeroclastic  solution.  The  addition  of  stores  did  not  affect  the  LCO  frequency,  which 
was  always  approximately  3.3  Hz.  In  addition  to  the  negative  offset,  under-wing 
stores  decreased  the  mode  1  amplitudes,  which  decreased  the  amount  of  bending. 
Tip  stores  increased  the  mode  1  amplitudes,  which  increased  the  amount  of  bending. 

Aerodynamic  store  forces  increase  LCO  amplitudes.  This  effect  was  demon¬ 
strated  when  the  aerodynamic  store  forces  were  magnified  in  order  to  better  discern 
their  effect  on  LCO.  Magnifying  the  forces  showed  that  the  twist  of  the  wing  was 
driven  by  mode  1  and  mode  2  coupling.  The  bending  force  primarily  affected  the 
LCO  by  directly  contributing  to  mode  1.  Magnifying  the  bending  force  increased 
the  amplitude  of  the  LCO.  This  was  expected  because  increasing  the  restoring  force 
resulted  in  greater  overshoot  beyond  the  neutral  position.  Magnifying  but  reversing 
the  direction  of  the  aerodynamic  store  forces  decreased  the  restoring  force  resulting 
in  less  overshoot  and  a  decreased  amplitude  LCO. 
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Interference  from  under-wing  stores  decreases  LCO  amplitudes.  The  larger  the 
diameter  of  the  under-wing  store,  the  larger  the  interference  with  the  flow  field  on  the 
bottom  of  the  wing.  This  led  to  increased  energy  quenching  and  a  lower  amplitude 
LCO.  Multiple  stores  were  added  to  the  Goland+  wing.  Again,  it  was  found  that 
the  larger  the  store,  the  more  energy  quenching  that  was  present  due  to  interference 
with  the  flow  field  on  the  bottom  of  the  wing. 

In  summary,  The  energy  quenching  mechanism  responsible  for  limit-cycle  os¬ 
cillation  in  the  Goland+  wing  has  been  determined.  The  primary  bending  and  tor¬ 
sional  modes  couple,  resulting  in  periodically  appearing  and  disappearing  lambda 
and  trailing-edge  shocks  which  decrease  the  restoring  forces  on  the  wing,  limiting 
the  flow  of  energy  from  the  fluid  into  the  structure.  Stores  affect  LCO  in  two  offset¬ 
ting  ways.  Adding  under-wing  stores  interferes  with  the  airflow  on  the  bottom  of  the 
wing  which  also  reduces  the  restoring  forces,  thereby  providing  additional  quenching 
of  the  LCO.  Under-wing  and  tip  stores  also  add  carriage  loads  into  the  wing  which 
increases  the  flow  of  energy  from  the  fluid  into  the  structure,  thereby  amplifying  the 
LCO. 

5.2  Recommendations  for  Future  Research 

The  Goland+  wing  was  analyzed  with  an  inviscid  fluid  solver  and  using  the 
thin-plate  spline  to  transfer  data  between  the  aerodynamic  and  structural  grids. 
These  splines  conserved  total  forces  and  moments  but  did  not  conserve  energy.  The 
next  step  in  the  research  should  involve  incorporation  of  an  energy-conserving  spline 
that  more  accurately  transfers  data  between  the  grids. 

FLUENT’s  grid  deformation  algorithms  were  the  main  limitation  preventing 
the  use  of  the  viscous  solver  and  FLLIENT’s  turbulent  models.  Grid  deformation 
could  be  controlled  externally  by  incorporating  a  faster  and  more  robust  grid  defor¬ 
mation  algorithm  tailored  for  deforming  wings.  This  could  be  implemented  through 
a  UDF,  bypassing  FLLIENT’s  internal  grid  deformation  algorithms.  FLLIENT  In- 
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corporated  is  also  addressing  this  short-coming  and  may  provide  a  faster  and  more 
robust  grid  deformation  algorithm  in  a  future  release.  With  a  faster,  more  robust 
grid  deformation  algorithm,  and  using  the  viscous  solver,  the  flow  around  the  wing 
tip  and  around  the  stores  could  be  examined  in  more  detail.  A  viscous  analysis 
into  how  stream-wise  movement  of  stores  could  then  be  conducted.  Viscous  analysis 
would  also  allow  for  nonlinearities  due  to  shock- wave/boundary-layer  interaction  to 
be  studied. 

Inertial  forces  were  added  to  the  Goland+  wing  in  order  to  compare  how  store 
mass  position  affected  LCO  as  opposed  to  the  store  aerodynamics.  In  this  research, 
inertial  forces  were  added  to  the  store  forces  by  simulating  the  addition  of  a  store 
mass.  The  structural  model  was  not  modified;  the  mass  matrix,  stiffness  matrix, 
and  the  modal  matrix  did  not  contain  any  store  mass  effects.  In  future  research, 
the  structural  model  should  be  modified  to  account  for  these  store  masses  in  order 
to  validate  the  results  what  were  obtained  by  simulating  the  masses.  The  effect  of 
adding  masses  at  different  span-wise  and  stream-wise  positions  on  LCO  and  on  the 
flutter  stability  boundary  should  be  further  studied. 

There  is  significant  interest  in  the  discovery  of  LCO  during  flight  test  of  an 
F-16.  Future  research  should  try  to  duplicate  this  flight  test  scenario.  Initial  simu¬ 
lations  of  this  scenario  did  not  result  in  LCO  (Appendix  F).  In  the  initial  inviscid 
simulations,  only  the  wing  and  the  tip  launcher  were  modelled.  The  fuselage  and 
under-wing  stores  were  not  modelled.  The  dominant  bending  and  twisting  modes 
did  not  couple  into  a  single  degree-of-freedom  flutter  mode  as  did  the  Goland+  wing. 
Additional  testing  should  be  conducted  to  determine  whether  LCO  can  be  obtained 
with  this  simplified  wing  or  whether  additional  grid  details,  such  as  the  fuselage  and 
under-wing  stores,  need  to  be  included.  Viscous  modelling  may  also  be  necessary  to 
duplicate  the  flight  test  results. 
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Appendix  A.  Spline  Routines 

For  aeroelastic  analysis,  there  are  two  grids;  an  aerodynamic  grid  on  the  surface  of 
the  body,  and  a  structural  grid  on  the  structural  members  of  the  body.  Data  needs  to 
be  transferred  between  these  two  grids,  but  they  normally  are  non-point  matching. 
In  this  case,  a  spline  matrix  can  be  used  to  transfer  data  between  the  grids. 

A.l  Thin- Plate  Spline 

The  thin-plate  spline  method  provides  a  means  to  characterize  an  irregular  sur¬ 
face  by  using  functions  that  minimize  fQ  \D2v\2,  a  functional  similar  to  the  bending 
energy  of  a  thin  plate  [18].  Duchon  [18]  proved  that  only  one  solution  exists  in  the 
form 

a{t)  =  ^2  ^ a\t  —  a]'2  ln  \t  —  °|  +  aot  +  /3  (A.l) 

aeA 

with  =  0  and  =  0. 

The  function  in  A.l  can  be  used  to  construct  a  function  S(x),  given  its  values 
at  a  set  of  N  discrete  “nodal”  values  [50] 

N 

S(x)  =  P  +  aQxx  +  aQyy  +  aQzz  +  ^  Xi\x  -  Xi\2  In \x  -  Xi\.  (A.2) 

i= 1 

The  coefficients  A;,  ft,  a0x ,  aQy,  and  aQz  are  determined  by  solution  of  the  minimiza¬ 
tion  problem  [50].  Applying  A.2  to  each  of  the  N  nodes  yields 

N 

Hi  =  P  +  aQxXi  +  a0yyi  +  a0zZi  +  ^  A 5 \xi  -  Xj |2  In  %  -  xj\ .  (A.3) 

3= 1 

Equation  A.3  is  subject  to  the  following  side  conditions: 

N  N  N  N 

y  )  A i  y  ^  A pX{  y  ^  Aj|/.j  y  ]  A iz^  O.  (A.4) 

i= 1  2=1  2=1  2=1 
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Writing  equations  A. 3  and  A. 4  in  matrix  notation  gives 


{<$}  =  [B]{A}  +  [i2]{a} 


The  matrix  [ R ]  is  defined  as 


The  column  matrices  {5},  {A},  and  {a}  are  defined  as 

Si 
S. 2 

Sn 


(A.5) 

(A.6) 

(A.7) 

(A.8) 

(A.9) 

(A. 10) 
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Ai 


(A}  = 


and 


w  = 


Atv 

Po 

aox 

a°y 


Solving  equations  A. 5  and  A. 6  for  {A}  and  {a}  gives 


(A. 11) 


(A. 12) 


{A}=  [B]-1  -  [S]-1^]  [[JR]T’[JB]-1[je]]-1  [JR]7'[B]-1  W  (A. 13) 


and 


M  =  [[JS]I'[S]-1[JS]]-1  [UnB]-1^}-  (A. 14) 

Letting  subscript  s  represent  the  structural  grid  and  subscript  a  represent  the  aero¬ 
dynamic  grid,  equation  A. 5  can  then  be  written  as 

{<54  =  [£U]{Aj  +  [i4]{aJ-  (A. 15) 

The  transformation  spline  matrix  [S]  is  needed  in  the  form 

R}  =  [S]{<5s}.  (A. 16) 


A-3 


Therefore,  {<5S}  can  be  factored  out  of  equation  A. 15  and  the  equations  can  be  solved 
for  [S], 


[A]  =  [[Bas]  [[5,]”1  -  [B^R.]  [[i2jr[5j"1[i2j]“1  [RsfiBs]-1' 

+[i4]  [[i2s]T[S8]-1[i?s]]"1  [i2a]T[Ss]"1]  .  (A. 17) 

The  elements  of  [Bas]  and  [Bs]  are  defined  as 

BaSij  =  r2aSijHr2aSij)  (A.  18) 

and 

=r»«ln(r»«)  (A-19) 

where 

r2aSij  =  {xai  ~  xSj )2  +  (yai  -  ySj  f  +  (zai  -  zSj)2  (A.20) 

and 

r%  =  ( xsi  -  xSj)2  +  (ySi  -  ySj)2  +  (zSi  -  zSj)2.  (A. 21) 

[R0]  and  [Rs]  are  defined  by  equation  A. 9  where  the  a  subscript  stands  for  the 
aerodynamic  grid  and  the  s  subscript  stands  for  the  structural  grid. 

Once  the  spline  transformation  matrix  [S]  is  generated,  displacements  and 
coordinates  of  the  displaced  aerodynamic  grid  can  be  computed  from  displacements 
of  the  structural  grid  with  the  following: 

N  =  [SIN,  (A.22) 

[q]  —  [q]o  +  [<5]  (A. 23) 
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where  [q\0  is  the  original  undeformed  grid.  The  grid  coordinate  matrices  are  defined 


Qx  i 

Qyi 

Qzi 

0.X2 

(hn 

Qz2 

Qxn 

%N 

Qzn 

(A. 24) 


In  addition,  the  transferal  of  forces  from  the  aerodynamic  grid  to  the  structural  grid 
can  be  achieved  by  the  transpose  of  [S']  [1], 


[fs]  = 

where  the  force  matrices  are  defined  as 


F 

1  xi 

F 

1  y  i 

FZ1 

F 

A  X2 

F 

1  V2 

fZ2 

F 

1  XN 

F 

1  Vn 

FZn 

(A. 25) 


(A. 26) 


A  limitation  of  the  thin-plate  spline  is  that  no  two  structural-grid  points  can 
be  located  at  the  same  coordinates,  nor  can  all  the  structural-grid  points  he  in  the 
same  plane  [1]. 

A. 2  Infinite-Plate  Spline 

The  infinite-plate  spline  method  [28]  is  based  upon  the  small  deflection  equation 
of  an  infinite  plate.  Mathematically,  it  is  a  two-dimensional  implementation  of  the 
thin-plate  spline  method.  The  changes  to  equations  A. 3,  A. 8,  A. 9,  A.  12,  A. 20,  and 
A. 21  are  as  follows: 

N 

Hi  =  (3  +  a0xXi  +  a0yy.i  +  ^  Xj \xt  -  xj\2  In  [%-%[,  (A. 27) 

3=1 
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4  =  (xi  -  xof  +  (Vi  -  Vj)2’ 


(A. 28) 


r2aSij  =  (xai  -  xSj)'2  +  (yai  -  ySj)%  (A.31) 

and 

r%  =  {xSi  ~  xSjf  +  (ySi  -  ySj)2.  (A. 32) 

The  remaining  equations  are  the  exact  same  as  formulated  for  the  thin-plate 
spline.  Equations  A. 24  and  A. 25  remain  the  same  since  the  spline  is  still  used  for  a 
three-dimensional  problem.  Since  the  aerodynamic  grid  is  three  dimensional,  when 
creating  the  [R0]  matrix  as  shown  in  equation  A. 29,  the  aerodynamic  grid  points  are 
projected  to  the  spline  plane  on  which  the  structural  grid  lies.  Or  in  other  words, 
the  ^-component  of  the  grid  coordinates  are  ignored. 

Unlike  the  thin-plate  spline,  the  structural-grid  points  can  lie  on  the  same  plane 
for  the  infinite-plate  spline.  However,  they  can  not  lie  in  a  line.  Like  the  thin-plate 
spline,  no  two  structural-grid  points  can  be  located  at  the  same  coordinates  [1]. 
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A. 3  Spline  Conclusions 

Any  spline  matrix  can  be  used  with  the  aeroelastic  program.  Two  have  been  im¬ 
plemented,  the  thin-plate  spline  and  the  infinite-plate  spline.  If  a  three-dimensional, 
structural  grid  is  used,  the  thin-plate  spline  provides  superior  results  [50].  If  the 
structural  grid  is  a  two-dimensional,  flat  plate,  the  infinite-plate  spline  must  be  used 
because  the  thin-plate  spline  is  invalid  if  all  the  points  lie  in  a  single  plane.  If  the 
structural  grid  is  based  on  a  beam  model  with  all  the  grid  points  in  a  straight  line, 
neither  of  these  two  methods  will  work  because  neither  method  is  valid  if  all  the 
points  lie  in  a  line. 
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Appendix  B.  Structural  Model 


The  semi-discrete  equations  of  motion  can  be  written  as 

[M]{i}  +  [C]{«}  +  [if]{U}  =  {F}  (B.l) 

where  [M\  is  the  mass  matrix,  [C]  is  the  damping  matrix,  [. K ]  is  the  stiffness  matrix, 
{F}  is  the  applied  forces,  {«}  is  the  modal  displacements,  {ft}  is  the  modal  velocities, 
and  {it}  is  the  modal  accelerations  [30]. 

One  of  the  most  widely  used  integration  methods  for  solving  this  initial-value 
problem  is  the  Newmark  method  [30],  which  consists  of  the  following  equations: 


([M]  +  7A t[C]  +  pAt2[K})  {itjn+i  =  {F}n+i  -  [C\  ({ft}„  +  (1  -  7)At{ii}n) 

f  At 2  \ 

-  [K]  (  {u}n  +  A t{u}n  +  — (1  -  2(3){u}n  J  ,  (B.2) 

{iljn+1  =  {u}n  +  At  [(1  -  7 ){u}n  +  7{fi}n+i]  ,  (B.3) 

At2 

{u}n+i  =  {u}n  +  At{u}n  +  —  [(1  -  2 (3){u}n  -  2p{u}n+1\ .  (B.4) 

If  7  =  },  this  implicit  method  is  second-order  accurate  [30].  This  method  is  stable 
if  At  <  ^  where  uj  is  the  maximum  natural  frequency  [30].  If  a  linear  acceleration 
is  assumed,  P  —  \  and  ft  =  3.464  [30]. 

Assuming  a  linear  acceleration  and  no  structural  damping,  equations  B.2,  B.3, 
and  B.4  simplify  to: 
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[M\  + 


{u}n+\  =  {F}n+ 1  -  [K] 


^{w}n  +  At{w}n  + 


(B.5) 


{u}n+ 1  =  {u}n  +  (({«}n  +  («}n+l)  ,  (B.6) 

and 

At2 

{u}n+ 1  =  {u}n  +  A t{u}n  +  —  (2 {u}n  +  {n}n+l)  •  (B.7) 

O 

Within  the  aeroclastic  routine,  forces  are  calculated  at  the  centroid  of  each 
aerodynamic  grid  wing  face.  These  forces  are  then  splined  to  the  structural  grid 
nodes.  A  generalized  force  vector  for  the  modal  system  is  then  calculated  from 

{jy  =  [t]T{F,}.  (B.8) 

{Fg}  is  the  generalized  force  and  moment  vector,  {Fs}  is  the  force  and  moment 
vector  for  the  structural  grid  nodes,  and  [<f>]  is  the  modal  matrix,  the  columns  of 
which  are  the  eigenvectors.  There  is  a  slight  phase  lag  between  the  fluid  solver  and 
the  structure.  The  forces  can  be  extrapolated  forward  in  time  to  try  and  mitigate 
this  phase  lag  if  desired.  The  extrapolated,  generalized  force  and  moment  vector  is 
calculated  as 

{ Fg }  n  =  2-5{F9}at  —  2.0{Fff}Ar_i  +  0.5{F9}tv-2-  (B.9) 

Once  the  forces  are  determined,  Equation  B.5  is  first  solved  for  the  generalized 
acceleration  vector  {w}n+1  using  an  LUP-decomposition  [11],  Equations  B.6  and 
B.7  are  then  solved  directly  to  get  the  generalized  velocity  vector  {u}n+\  and  the 
generalized  displacement  vector  {;u}n+1 .  Displacement  and  rotation  of  the  structural 
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grid  nodes  are  then  found  from 


{4}  =  [*]{«}■*-!• 


(B.10) 


These  new  displacements  are  then  splined  to  get  new  aerodynamic  grid  dis¬ 
placements  and  coordinates  in  order  to  deform  the  mesh. 
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Appendix  C.  Grid  Convergence  Study 

Once  the  aeroelastic  program  was  validated,  it  was  used  to  study  LCO  of  the  Goland+ 
wing  described  in  section  3.2.2.  Four  grids  were  built  to  look  at  grid  convergence; 
a  coarse  grid  consisting  of  68,949  tetrahedral  cells  with  518  cell  faces  on  the  wing, 
a  medium  grid  consisting  of  269,596  tetrahedral  cells  with  13,358  cell  faces  on  the 
wing,  a  fine  grid  consisting  of  660,347  tetrahedral  cells  with  42,102  cell  faces  on  the 
wing,  and  a  viscous  grid  consisting  of  2,203,065  tetrahedral,  prism,  and  pyramid  cells 
with  54,886  cell  faces  on  the  wing.  The  coarse,  medium,  and  fine  grids  were  tested 
using  the  inviscid  solver.  The  viscous  grid  was  tested  using  the  Spalart-Allmaras 
turbulence  model.  These  grids  were  used  to  compute  steady  state  solutions  for 
Mach  0.92,  600  ft/sec  (182.88  m/sec),  at  angles  of  attack  of  0  to  20  degrees  in  2  degree 
increments.  Cl  and  Cm  plots  are  shown  in  Figure  C.l.  Below  11  degrees  angle-of- 
attack,  all  of  the  inviscid  grids  produced  similar  Cl  and  Cm  results.  The  viscous 
grid  Cl  results  are  less  and  the  Cm  results  are  more  than  the  results  produced  by  the 
inviscid  grids.  Above  11  degrees,  the  coarse  grid  is  under- predicting  CL  and  slightly 
over  predicting  Cm  while  the  viscous  results  and  medium  grid  inviscid  results  were 
approximately  the  same.  The  viscous  results  did  not  predict  a  stall  under  20  degrees, 
however,  a  small  separation  bubble  behind  the  leading  edge  did  appear  at  10  degrees 
angle-of-attack  and  continued  to  grow  larger  as  the  angle-of-attack  was  increased. 
The  inviscid  fluid  solver  failed  in  this  region  for  the  fine  grid  at  angles-of-attack 
greater  than  12  degrees. 

A  slice  through  the  wing  was  taken  at  3.2808  ft  (1  m)  span  and  the  coefficient 
of  pressure  for  the  airfoil  was  plotted.  As  can  be  seen  in  figures  C.2  -  C.12,  the 
coarse  grid  did  not  accurately  capture  the  pressure  distribution  across  the  wing.  At 
zero  degrees  angle-of-attack,  the  shock  was  located  at  a  chord  position  of  4.921  ft 
(1.5  m).  For  the  inviscid  grids,  as  the  angle-of-attack  was  increased  to  2  degrees, 
the  shock  moved  immediately  to  the  trailing  edge  of  the  wing.  For  the  viscous  grid, 
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C,  vs  Angle  of  Attack 
Goland+  Wing,  Mach  0.92,  600  ft/sec 


(a)  CL 


C,,  vs  Angle  of  Attack 
Goland+  Wing,  Mach  0.92,  600  ft/sec 


(b)  CM 

Figure  C.l  Cl  and  Cm  vs  Angle-of-Attack 
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0.0  Degrees  Angle  of  Attack 


Figure  C.2  CP  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  0  degrees 


2.0  Degrees  Angle  of  Attack 


Figure  C.3  Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  2  degrees 
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4.0  Degrees  Angle  of  Attack 


Figure  C.4  CP  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  4  degrees 


6.0  Degrees  Angle  of  Attack 


Figure  C.5  Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  6  degrees 
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8.0  Degrees  Angle  of  Attack 


Figure  C.6  CP  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  8  degrees 


1 0.0  Degrees  Angle  of  Attack 


Figure  C.7  Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  10  degrees 
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1 2.0  Degrees  Angle  of  Attack 


Figure  C.8  CP  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  12  degrees 


1 4.0  Degrees  Angle  of  Attack 


Figure  C.9  Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  14  degrees 
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Figure  C.10 


Figure  C.ll 


1 6.0  Degrees  Angle  of  Attack 


CP  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  16  degrees 


1 8.0  Degrees  Angle  of  Attack 


Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  18  degrees 
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20.0  Degrees  Angle  of  Attack 


Figure  C.12  Cp  plots  for  Mach  0.92,  600  ft/sec  (182.88  m/sec),  20  degrees 

the  shock  slowly  moved  aft  and  did  not  reach  the  trailing  edge  until  approximately 
8  degrees.  At  14  degrees  angle-of-attack  and  above,  the  viscous  results  showed  that 
the  flow  separated  behind  the  shock  on  the  top  of  the  wing. 

The  coarse  grid  failed  to  accurately  capture  the  pressure  distribution  at  the 
leading  and  trailing  edge  of  the  wing  but  it  did  capture  the  pressure  distribution  in 
the  center  of  the  wing.  A  fourth  inviscid  grid  was  then  built  that  fell  between  the 
coarse  and  medium  grids.  This  grid  correctly  captured  the  pressure  distribution  at 
the  leading  and  trailing  edge  of  the  wing,  but  with  fewer  cells  than  the  medium  grid. 
It  consisted  of  194,780  tetrahedral  cells  with  9,178  cell  faces  on  the  wing.  This  grid 
was  used  to  compute  a  clean-wing  baseline  solution  against  which  all  other  LCO 
simulations  were  compared. 

At  600  ft/sec  (182.88  m/sec)  and  Mach  0.92,  the  dynamic  pressure  is  beyond 
the  flutter  point  and  was  shown  by  Snyder,  et  al.  [51]  to  result  in  LCO.  For  the 
Goland+,  these  flow  conditions  fall  within  the  transonic  flutter  dip  region.  This 
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case  was  used  as  a  baseline  for  all  simulations  in  this  study.  To  match  the  condi¬ 
tions  of  Snyder,  et  al.  [51],  the  far-held  density,  temperature,  and  pressure  were  set 
to  0.0023771  slugs/ ft3  (1.225  kg/m3),  518.67° R  (288.15  K ),  and  722.1813  lb/ ft2 
(34578.04  Pa)  respectively,  while  R  was  set  to  585.7438  ( R  =  97.951  ) . 

LCO  for  the  Goland+,  clean- wing  was  computed  on  all  five  grids  discussed 
above,  at  600  ft/sec  (182.88  m/sec),  Mach  0.92.  The  three  coarsest  inviscid  grids 
correctly  captured  the  LCO.  This  was  initially  surprising,  because  the  coarse  grid 
did  not  accurately  capture  the  pressure  distribution  on  the  wing  at  the  leading  and 
trailing  edges,  nor  did  it  accurately  define  the  shock  location  due  to  the  large  size 
of  the  cells.  However,  because  the  structural  grid  was  coarser  than  even  the  coarse 
aerodynamic  grid,  all  of  the  aerodynamic  grids  produced  a  similar  force  distribu¬ 
tion  when  transferred  to  the  structural  grid.  The  amplitudes  of  the  LCO  and  the 
modal  contributions  were  approximately  the  same  regardless  of  the  grid.  The  only 
significant  variation  in  the  LCO  due  to  grid  refinement  was  in  the  frequency.  The 
frequency  was  3.37  Hz  for  the  coarse  grid,  3.31  Hz  for  the  baseline  grid,  and  3.27  Hz 
for  the  medium  grid.  The  computation  time  varied  greatly.  It  took  approximately 
5.5  hours  to  produce  one  second  of  data  using  six  2.2  GHz  opteron  processors  on  the 
coarse  grid,  15.6  hours  on  the  baseline  grid,  and  21  hours  on  the  medium  grid. 

The  fine  grid  produced  only  6.584  seconds  of  data  before  the  fluid  solver  failed. 
This  failure  was  due  to  the  twist  producing  high  angles-of-attack  at  the  wing  tip,  and 
the  solver’s  inability  to  handle  the  adverse  pressure  gradient  that  appeared  behind 
the  leading  edge.  This  failure  was  in  the  same  location  that  the  steady,  viscous 
runs  showed  a  very  small  separation  bubble.  The  coarse  grids  averaged  out  this 
small  problem  area  so  that  the  solver  did  not  crash.  When  the  solver  failed,  Cl  was 
oscillating  at  a  frequency  of  3.15  Hz  and  the  frequency  was  increasing.  With  this 
grid,  it  took  27.3  hours  to  produce  one  second  of  data  using  six  2.2  GHz  opteron 
processors. 
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The  viscous  grid  using  the  Spalart-Allmaras  turbulence  model  was  run  for  less 
than  3  seconds  because  of  very  large  computation  times.  With  this  grid,  it  was 
taking  seven  days  to  produce  one  second  of  data  using  sixteen  2.2  GHz  opteron  pro¬ 
cessors.  As  the  amplitudes  continued  to  grow  and  the  motion  continued  to  increase, 
the  computational  time  required  for  each  time  step  was  increasing.  Therefore,  this 
aeroclastic  analysis  tool  was  found  to  be  inadequate  for  viscous  calculations. 
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Appendix  D.  GolandA  with  NACA  65-004  Airfoil 

To  determine  whether  the  airfoil  affects  the  flow  of  energy  between  the  fluid  and  the 
structure,  or  whether  the  flow  of  energy  was  independent  of  the  airfoil,  the  Goland+, 
biconvex  airfoil  was  replaced  with  a  NACA  65-004  airfoil.  Since  the  airfoil  affects 
the  shock  strength,  location,  and  movement,  it  was  expected  to  affect  the  amplitude 
of  the  LCO.  With  the  biconvex  airfoil,  LCO  was  established  by  6.4  seconds  with  a 
mode  1  amplitude  of  ±8.8  as  discussed  above.  With  the  same  initial  perturbation 
and  flow  conditions,  the  NACA  65-004  airfoil  mode  1  amplitude  only  realized  ±1.4 
after  15  seconds  (Figure  D.l).  This  demonstrated  that  the  airfoil  was  critical  to  the 
transfer  of  energy  between  the  structure  and  fluid.  A  large  initial  velocity  pertur¬ 
bation  was  then  given  to  the  NACA  65-004  configuration  (Figure  D.2).  The  large 
perturbation  quickly  accelerated  the  motion  past  the  LCO  conditions  seen  with  the 
bi-convex  wing.  Mode  1  and  mode  2  were  still  in  phase  as  they  were  for  the  bi¬ 
convex  wing.  This  led  to  the  same  trailing-edge  and  leading-edge  shock  structure 
(Figure  D.3)  as  previously  seen.  At  6.365  seconds  the  inviscid  fluid  solver  crashed 
due  to  the  large  anglc-of-attack.  Based  on  the  shock  structure  and  mode  coupling, 
it  was  hypothesized  that  it  would  settle  into  LCO,  although  at  an  angle-of-attack 
higher  than  the  inviscid  fluid  solver  was  capable  of  simulating.  Regardless,  this  case 
demonstrated  the  importance  of  the  airfoil  on  the  LCO  solution.  The  airfoil  dic¬ 
tated  how  fast  the  oscillations  grew  and  the  final  LCO  amplitudes.  The  coupling 
of  the  modes  and  the  corresponding  shock  motion  determined  whether  a  damping 
mechanism  was  present  that  could  lead  to  LCO. 
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Figure  D.2 


NACA  65-004 


Figure  D.l  NACA  65-004  Cl  and  Modal  Amplitudes 


NACA  65-004  Airfoil 
Mach  0.92,  600  ft/sec 


NACA  65-004  (Large  Velocity  Perturbation)  Cl  and  Modal  Ampli¬ 
tudes 
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Figure  D.3 
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Appendix  E.  LCO  at  Positive  Angle- of- Attack 

To  determine  whether  the  angle-of-attack  of  the  wing  affects  the  flow  of  energy  be¬ 
tween  the  fluid  and  the  structure,  the  root  angle-of-attack  was  changed  from  0.0  de¬ 
grees  to  2.0  degrees.  Changing  the  angle-of-attack  changed  the  static  aeroelastic 
solution,  causing  a  displacement  of  the  static  values  about  which  the  modes  and 
Cl  (Figure  E.l)  oscillated.  The  wing  oscillated  about  a  Cl  of  0.1676  instead  of  0.0 
and  about  a  tip  angle-of-attack  of  1.36  degrees  instead  of  0.0  degrees.  This  new 
neutral  point  was  the  new  static  aeroelastic  solution  for  the  wing  at  2.0  degrees 
angle-of-attack.  The  pattern  of  the  energy  flow  between  the  fluid  and  structure  was 
unaffected. 
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Goland+  Clean  Wing 
Mach  0.92,  600  ft/sec 
2  degrees  Angle  of  Attack 


Time 


Figure  E.l  Cl  and  Modal  Amplitudes  at  2.0  Degrees  Angle-of- Attack 
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Appendix  F.  F-16  Results 

Air  Force  SEEK  EAGLE  Office  used  a  block  40  F-16C  for  LCO  testing  [15].  The 
F-16  has  a  cropped  delta  wing,  blended  with  the  fuselage.  It  has  a  wingspan  of 
32  ft,  8  in  (9.9568  m)  and  a  NACA  64A-204  airfoil.  The  wing  aspect  ratio  is  3.2  and 
the  leading-edge  sweep  is  40  degrees  [15].  The  configuration  shown  in  Figure  F.l 
exhibited  LCO  in  level  flight  test  conditions  ranging  from  Mach  0.70  to  Mach  0.95 
at  2,000  ft,  5,000  ft,  and  10,000  ft  pressure  altitude  [15]. 

An  attempt  to  compute  the  F-16  LCO  state  with  the  aeroelastic  program  was 
made.  The  store  configuration  and  mass  properties  used  by  the  modal  structural 
model  obtained  from  Denegri  [15],  are  shown  in  Tables  F.l  and  F.2.  The  modal 
structural  model  consisted  of  26  modes  which  ranged  in  frequency  from  5.21  Hz  to 
24.33  Hz  [15].  The  first  six  modes  are  shown  in  Figure  F.2.  An  inviscid  grid  made  up 
of  1,189,805  tetrahedral  cells  was  created  for  this  problem.  The  under-wing  stores 
and  the  fuselage  were  not  modelled  in  the  aerodynamic  grid,  but  the  tip  launcher 
was  modelled.  The  wing  was  projected  to  the  centerline  through  the  area  located  by 
the  fuselage,  therefore,  the  total  span  was  correct  (Figure  F.3).  Denegri  [15]  found 
that  the  dominant  LCO  modes  were  anti-symmetric,  therefore,  both  sides  of  the 
F-16  wing  were  modelled.  The  2,000  ft  pressure  altitude,  Mach  0.90  test  condition 
was  simulated.  Acceleration  data  for  the  aerodynamic  grid  point  corresponding  to 
the  forward,  vertical  accelerometer  was  collected.  This  was  compared  to  the  ffight- 


Figure  F.l  F-16C  Store  Configuration 
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Table  F.l  F-16  Store  Configuration  and  Attachment  Reference  Points 


Station 

Location 

x  (in.) 

y  (in.) 

z  (in.) 

Store 

Suspension  equipment 

1 

Wingtip 

380.46 

180.00 

92.00 

None 

LAU-129/A  launcher 

2 

Underwing 

371.56 

157.00 

91.35 

AIM-9P  missile 

LAU-129/A  launcher 

3 

Underwing 

349.67 

120.00 

90.72 

Air-Ground  missile 

Launcher  /  pylon 

4 

Underwing 

316.87 

71.00 

90.02 

Empty  370  gal.  tank 

Fuel  Pylon 

Left  wing  station  loading.  For  the  flight  test,  the  right  wing  was  configured  as  a  mirror  image 
of  the  left  wing. 


Table  F.2  F-16  Store  Mass  Properties 


Store 

Weight 

Center  of  gravity"1 

Moments  of  inertia 
(slug-ft.2) 

(lb.) 

x(in.) 

y(in.) 

z(in.) 

Roll 

Pitch 

Yaw 

LAU-129/A  wingtip  launcher 

84.9 

-15.20 

0.00 

-0,01 

- 

13.73 

13.70 

AIM-9P  missile,  LAU-129/A  un¬ 
derwing  launcher,  and  missile  py¬ 
lon 

276.3 

-14.65 

0.01 

-13.30 

1.31 

68.63 

67.77 

Air-ground  missile,  launcher,  and 
weapon  pylon 

898.7 

-2.27 

0.02 

-20.72 

19.72 

134.74 

116.95 

370  gallon  fuel  tank  (empty)  and 
pylon 

470.5 

-1.85 

-0.18 

-10.00 

197.93 

187.64 

Positive  values  are  aft,  outboard,  or  above  the  store  attachment  reference  point  x-y-z  locations 
(see  Table  F.l). 


test  results  shown  in  Figure  F.4  [15].  This  figure  shows  that  the  wing-tip  response 
amplitude  was  about  ±3g  or  1, 158  in/sec2  (29.42  m/s2)  for  this  test  condition. 

With  an  initial  modal  displacement  perturbation  of  0.02  to  mode  2,  after  11  sec¬ 
onds,  the  wing-tip  response  amplitude  was  around  180  in/sec 2  (4.57  m/s2),  but  still 
growing.  Since  the  computation  time  for  this  grid  was  over  50  hours  for  one  second 
of  data  on  sixteen  2.2  GHz  Opteron  processors,  it  would  take  too  long  to  determine 
whether  LCO  would  develop  using  this  small  initial  perturbation.  The  accelerations 
had  a  frequency  of  8.4  Hz,  which  was  close  to  that  reported  by  Denegri  [15]. 

A  large,  initial-velocity  perturbation  of  10.0  was  given  to  mode  4  to  determine 
if  run  times  could  be  reduced.  The  large  initial  perturbation  kicked  the  wing  to 
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(a)  Mode  1  (5.212726  Hz)  (b)  Mode  2  (8.167878  Hz)  (c)  Mode  3  (8.466704  Hz) 


(d)  Mode  4  (8.671919  Hz)  (e)  Mode  5  (10.56468  Hz)  (f)  Mode  6  (10.89219  Hz) 

Figure  F.2  Mode  Shapes  for  F-16C 


Figure  F.3  F-16C  Inviscid  Grid 
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Figure  F.4  F-16C  Wing-Tip  Response 

the  acceleration  levels  seen  in  flight  test,  but  did  not  result  in  LCO.  The  amplitudes 
continued  to  grow  resulting  in  divergent  flutter  (Figure  F.5).  Modes  2  and  4  were  the 
dominant  modes.  Mode  2  was  the  first  anti-symmetric  bending  mode  and  mode  4 
was  the  first  anti-symmetric  twisting  mode.  These  two  modes  were  107  degrees 
out  of  phase  and  did  not  couple  into  a  single  degree-of-freedom  flutter  mode  as  did 
the  Goland+  wing  (Figure  F.6).  These  results  imply  that  the  structural  model,  as 
implemented,  did  not  correctly  or  completely  model  the  structure  or  that  there  is  a 
different  mechanism  responsible  for  LCO  in  the  F-16  than  that  found  in  the  Goland+ 
wing.  The  different  mechanism  could  be  boundary-layer/shock-wave  interaction  or 
shock-induced  separation  on  the  wing,  both  of  which  require  viscous  modelling. 

Additional  testing  is  required  to  determine  whether  LCO  can  be  obtained  with 
this  simplified  wing  or  whether  additional  grid  details,  such  as  the  fuselage  and 
underwing  stores,  need  to  be  included.  Viscous  modelling  may  also  be  necessary  to 
duplicate  the  flight  test  results.  Other  researchers  are  also  trying  to  simulate  this 
flight-test  scenario  [14,  54,  55],  but  to  date,  no  one  has  been  able  to  completely 
duplicate  the  flight-test  results. 
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Figure  F.5  F-16  Wing-Tip  Acceleration  (Large  Initial  Condition) 
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